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A review of the treatment of boundaries in general relativity is presented with the emphasis on 
application to the formulations of Einstein's equations used in numerical relativity. At present, it 
is known how to treat boundaries in the harmonic formulation of Einstein's equations and a tetrad 
formulation of the Einstein-Bianchi system. However, a universal approach valid for other formula- 
tions is not in hand. In particular, there is no satisfactory boundary theory for the 3+1 formulations 
fSJ ' which have been highly successful in binary black hole simulation. I discuss the underlying problems 

T-H , that make the initial-boundary value problem much more complicated than the Cauchy problem. I 

c7^ ' review the progress that has been made and the important open questions that remain. 
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^ \ PACS numbers: PACS number(s): 04.20.-q, 04.20.Cv, 04.20.Ex, 04.25.D- 

^^ ' Science is a differential equation. Religion is a boundary condition. (Alan Turing, quoted in J.D. Barrow, "Theories 

t:]" . of Everything" ) . 
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"! I. INTRODUCTION 

^ ■ There are no natural boundaries for the gravitational field analogous to the conducting boundaries that play a major 

OX)| role in electromagnetism. In principle, there is no need to introduce any. The behavior of the universe as a whole can 

be posed as an initial value (Cauchy) problem. In an initial value problem, data is given on a spacelike hypersurface 

rvq ' Sq. The problem is to determine a solution in the future domain of dependence 'D^{So), which consists of those 

^ , points whose past directed characteristics all intersect Sq. The problem is well-posed if there exists a unique solution 
TJ" ' which depends continuously on the initial data. The pioneering work of Y. Bruhat |l|| showed that the initial value 
'^ ' problem for the (classical) vacuum gravitational field is well-posed. Assuming that matter fields do not spoil things, 
this suggests that the global cosmological problem of treating the universe as a whole can be solved in a physically 
meaningful way, i.e. in a way such that the solution does not undergo uncontrolled variation under a perturbation 
CO : of the initial data. This is indeed the case for the presently accepted cosmological model of an accelerating universe 
^^ (positive cosmological constant) where the conformal boundary at future null infinity Z+ is spacelike. In a conformally 

^^ , compactified picture, I'^ acts as a spacelike cap on the future evolution domain and no boundary condition is necessary 

. . ' or indeed allowed. 

K*" ', In practice, of course, treating an isolated system as part of a global cosmological spacetime is too complicated 

k>( • a problem without oversimplifying assumptions such as isotropy or homogeneity. One global approach applicable 
''»_] \ to isolated systems is to base the Cauchy problem on the analogue of a foliation of Minkowski spacetime by the 

d • hyperboloidal hypersurfaces 

f -x'^ -y^ - z'^ = T'^, t>T. (1.1) 

In a Penrose conformally compactified picture [2], this foliation asymptotes to the light cones and extends to a 
foliation of future null infinity I^. The analogue in curved spacetime is a foliation by positive constant mean 
curvature hypersurfaces. Since no light rays can enter an asymptotically flat spacetime through X"*", no boundary 
data are needed to evolve the interior spacetime. In addition, the waveform and polarization of the outgoing radiation 
can be unambiguously calculated at T"*" in terms of the Bondi news function 'S'l . This approach was first extensively 
developed by Friedrich 4] who formulated a hyperbolic version of the Einstein-Bianchi system of equations, which 
is manifestly regular at Z+, in terms of the conformally rescaled metric, connection and Weyl curvature. This is 
potentially the basis for a very attractive numerical approach to simulate globalproblems such as gravitational wave 
production. For reviews of progress on the numerical implementation see [5|-l7|. There has been some success in 
simulating model axisymmetric problems ,8]. More recently, there have been other attempts at the hyperboloidal 
approach based upon the Einstein equations for the conformal metric. Zenginoglu [9| has implemented a code based 
upon a generalized harmonic formulation in which the gauge source terms produce a hyperbolic foliation. A mixed 
hyperbolic-elliptic system proposed by Moncrief and Rinne [10| has been implemented as an axisymmetric code [11| 
which produces long term stable evolutions. Another hyperbolic-elliptic system based upon a tetrad approach has 
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been developed by Bardeen, Sarbach and Buchman [1^ . In spite of the attractiveness of the hyperboloidal approach 
and its success with model problems, considerable work remains to make it applicable to systems of astrophysical 
interest. 

A different global approach is to match the Cauchy evolution inside a a finite worldtube to an exterior characteristic 
evolution extending to T^ . In this approach, called Cauchy-characteristic matching, the characteristic evolution is 
constructed by using the Cauchy evolution to supply characteristic data on an inner worldtube, while the characteristic 
evolution supplies the outer boundary data for the Cauchy evolution. The success of Cauchy-characteristic matching 
depends upon the proper mathematical and computational treatment of the initial-boundary value problem (IBVP) 
for the Cauchy evolution. This approach has been successfully implemented in the linearized regime [1J| but also 
needs considerable additional work to apply to astrophysical systems. See [lj| for a review. 

Instead of a global treatment, the standard approach in numerical relativity, as in computational studies of other 
hyperbolic systems, is to introduce an artificial outer boundary. Ideally, the outer boundary treatment is designed to 
represent a passive external universe by allowing radiation to cross only in the outgoing direction. This is the primary 
application of the IBVP. Other possible applications, which I will not consider, are the timelike conformal boundary 
to a universe with negative cosmological constant and the membranes which play a role in higher dimensional theories. 
While there are no natural boundaries in classical gravitational theory, boundaries do play a central role in the ideas of 
holographic duality introduced in higher dimensional attempts at quantum gravity. Such applications are also beyond 
the scope of this review as well as beyond my own expertise. Here, I confine my attention to 4-dimensional spacetime, 
although the techniques governing a well-posed IBVP readily extend to hyperbolic systems in any dimension. 

In the IBVP, data on a timelike boundary T, which meets So in a surface So, are used to further extend the solution 
of the Cauchy problem to the domain of dependence V'^ {Sq U T) . In the simulation of an isolated astrophysical system 
containing neutron stars and black holes, the outer boundary T is coincident with the boundary of the computational 
grid and Sq is topologically a sphere surrounding the system. However, for purposes of treating the underlying 
mathematical and computational problems, it suffices to concentrate on the local problem in the neighborhood of 
some point on the intersection Bq between the Cauchy hypersurface Sq and the boundary T. For hyperbolic systems, 
the global solution in the spacetime manifold M can be obtained by patching together local solutions. This is because 
the finite speed of propagation allows localization of the problem. The setting for this local problem is depicted in 

Fig.m 




FIG. 1: Data on the S-manifolds So and T, which intersect in the 2-surface Bo, locally determine a solution in the spacetime 
manifold JVi 



The IBVP for Einstein's equations only received widespread attention after its importance to the artificial outer 
boundaries used in numerical relativity was pointed out [l5| . The first strongly well-posed IBVP was achieved for a 
tetrad version of the Einstein-Bianchi system, expressed in first differential order form, which included the tetrad, 
connection and curvature tensor as evolution fields [16| . Strong well-posedness guarantees the existence of a unique 
solution which depends continuously on both the Cauchy data and the boundary data. A strongly well-posed IBVP 
was later established for the harmonic formulation of Einstein's equations as a system of second order quasilinear 
wave equations for the metric [13] • The results were further generalized in [1^ [I^ to apply to a general quasilinear 
class of symmetric hyperbolic systems whose boundary conditions have a certain hierarchical form. 



A review of the IBVP in general relativity must per force be of a different nature than a review of the Ca uchy 
problem. The local properties of the Cauchy problem are now well understood. Several excellent reviews exist |20l - [23| . 
For the IBVP, the results are not comprehensive and are closely tied to the choice of hyperbolic reduction of Einstein's 
equations. There are only a few universal features and, in particular, there is no satisfactory treatment of the 3+1 
formulation which is extensively used in numerical relativity. For that reason, I will adopt a presentation which differs 
somewhat from the standard approach with the motivation of setting up a bare bones framework whose flexibility 
might be helpful in further investigations. 

My presentation is also biased by the important role of the IBVP in numerical relativity, which treats Einstein's 
equations as a set of partial differential equations (PDEs) governing the metric in some preferred coordinate system. 
On the other hand, from a geometrical perspective, one of the most fundamental and beautiful results of general 
relativity is that the properties of the local Cauchy problem can be summed up in geometric terms independent of 
any coordinates or explicit PDEs. This geometric formulation only came about after the Cauchy problem was well 
understood from the PDF point-of-view. The importance of the geometric approach to the numerical relativist is 
that it supplies a common starting point for discussing and comparing different formulations of Einstein's equations. 
Presently, the PDF aspects of metric formulations of the IBVP are only understood in the harmonic formulation. In 
order to transfer this insight into other formulations a geometric framework can serve as an important guide. For 
that reason, I will shift often between the PDF and geometric approach. When the emphasis is on the geometric side 
I will use abstract indices, e.g. w" to denote a vector field, and on the PDF side I will use coordinate indices, e.g. 
yfi _ (^^*^^«)^ to denote components with respect to spacetime coordinates a;^ = (^,a;*). 

The standard mathematical approach to the IBVP is to first establish the well-posedness of the underlying Cauchy 
problem, and next the local half-space problem. If these individual problems are well-posed (in a sense to be qualified 
later) then the problems with more general boundaries will also be well-posed. Thus I start my review in Sec. Ullbv 
first providing some brief background material for the Cauchy problem. 

Next, in Sec. IIIIl I point out the complications in going from the Cauchy problem to the IBVP. The IBVP for 
Einstein's equations is not well understood due to problems arising from the constraint equations. The motivation 
for this work stems from the need for an improved understanding and implementation of boundary conditions in the 
computational codes being used to simulate binary black holes. The ability to compute the details of the gravitational 
radiation produced by compact astrophysical sources, such as coalescing black holes, is of major importance to the 
success of gravitational wave astronomy. If the simulation of such systems is based upon a well-posed Cauchy problem 
but not a well-posed IBVP then the results cannot be trusted in the domain of dependence of the outer boundary. In 
See's llVi IVl and IVll I present the underlying mathematical theory. 

Early computational work in general relativity focused on the Cauchy problem and the IBVP only received consid- 
erable attention after its importance to stable and accurate simulations was recognized. I discuss some history of the 
work on the IBVP in Sec. IVIIi for the purpose of pointing out some of the partial successes and ideas which may be 
of future use. 

In addition to the mathematical issue of an appropriate boundary condition, the description of a binary black hole 
as an isolated system raises the physical issue of the appropriate outer boundary data. In the absence of an exterior 
solution, which could provide this data by matching, the standard practice is to set this data to zero. This raises 
the question, discussed in Sec. I Villi of how to formulate a non-refiecting outer boundary condition in order to avoid 
spurious incoming radiation. 

See's IIXI and 1x1 describe the two strongly well-posed formulations of the IBVP, which are known at the present time. 
Neither is based upon a 3 + 1 formulation and they both handle the constraints in different ways. This prompts the 
discussion, in Sec. lXIl of constraint enforcement in the 3-1-1 formulations. The resolution of issues regarding geometric 
uniqueness, discussed in Sec. IXIIl would shed light on the universal features of the IBVP that would perhaps guide 
the way to a successful 3-1-1 treatment. 

There is also the computational problem of turning a well-posed IBVP into a stable and accurate evolution code. 
I will not go into the details of the large range of techniques which are necessary for the successful implementation 
of a numerical relativity code. Since initiating this review, I have learned of a separate review in progress [24!| which 
covers such numerical techniques in great detail. Building a numerical relativity code is a complex undertaking. 
As observed by Post and Votta [25[ in a study of large scale computational projects, "the peer review process in 
computational science generally doesn't provide as effective a filter as it does for experiment or theory. Many things 
that a referee cannot detect could be wrong with a computational science paper. . . The few existing studies of error 
levels in scientific computer codes indicate that the defect rate is about seven faults per 1000 lines of Fortran" . They 
emphasize that "New methods of verifying and validating complex codes are mandatory if computational science is to 
fulfill its promise for science and society" . These observations are especially pertinent for numerical relativity where 
validation by agreement with experiment is not yet possible. In that spirit, I discuss the code tests that have been 
proposed and carried out for the gravitational IBVP in Sec. IXIIIl 

My aim has been to present the background material which might open new avenues for a better understanding of 



the IBVP and lead to progress on some of the important open questions posed in Sec. IXIVI 

II. THE CAUCHY PROBLEM 

Here I summarize those aspects of the Cauchy problem which are fundamental to the IBVP. For more detail, 
see [13-111. 

In contrast to Newtonian theory, which describes gravity in terms of an elliptic Poisson equation that propagates 
the gravitational field instantaneously, the retarded interactions implicit in general relativity give rise to new features 
such as gravitational waves. Wave propagation results from the mathematical property that Einstein's equations can 
be reduced to a hyperbolic system of PDEs. However, the coordinate freedom in Einstein's theory admits gauge waves 
which propagate with arbitrarily high speeds, including speeds faster than light. Einstein's equations are not a priori 
a hyperbolic system in which propagation speeds must be bounded and for which an initial value problem can be 
posed. 

This crucial step in going from Einstein's equations to a hyperbolic system has been highlighted by Friedrich as 



the process of hyperbolic reduction 26[. The first and most famous example of hyperbolic reduction was through 
the introduction of harmonic coordinates, which led to the classic result that the Cauchy problem for the harmonic 
formulation of Einstein's equations is well-posed [l|. Here I summarize the hyperbolic reduction of Einstein's equations 
in in terms of generalized harmonic coordinates x" = (t, a;') = (i, x, y, z), which are functionally independent solutions 
of the curved space scalar wave equation 

□x" = -La„(y^g"^9^x'^) = -f^ (2.1) 

where F^ are gauge source functions 26] . In terms of the connection F^ „ , these harmonic conditions are 

C^:=F^-f^ = 0, (2.2) 



where 



r^ = a'^^Kp = --^d^i./^gg'^^). (2.3) 

The hyperbolic reduction of the Einstein tensor results from setting 

E^"" := C" - W^T'> + -g'^'^S/pC = 0, (2.4) 

where C" is treated formally as a vector field in constructing the "covariant" derivatives V^C. (In generalized 
harmonic formulations based upon a background connection, C^ is a legitimate vector field. See Sec. 1X1) 

When the harmonic gauge source functions have the functional dependence T'^{x,g), the principal part of (j2.4p 
reduces to the wave operator acting on the densitized metric, i.e. 

E^"' = -^=da \g°'^dii{yf^g^"')\ + lower order terms. (2.5) 

Thus the harmonic evolution equations (j2.4|) are quasilinear wave equations for the components of the densitized 
metric \/—gg^'^ ■ The well-posedness of the Cauchy problem for the system (J2.4I) then follows from known results for 
systems of quasilinear wave equations. (It is important to bear in mind that such results are local in time since there 
is no general theory for the global existence of solutions to nonlinear equations.) 

In turn, the well-posedness of the Cauchy problem for the harmonic Einstein equations also follows provided that 
the harmonic conditions C^ = are preserved under the evolution. The proof of constraint preservation results from 
applying the contracted Bianchi identity V^G'^'" = to (12. 4p . This leads to a homogeneous wave equation for C^, 

VVpC^ + i^^C^O. (2.6) 

If the initial data enforce 

nso^O (2.7) 



and 

dtC^\so = (2.8) 



then the unique solution of p.6p is C = 0. It is easy to satisfy ()2.7|) by algebraicahy determining the initial values of 
dtg^* in terms of the initial values of g^^^ and their spatial derivatives. In order to see how to satisfy (|2.8|) note that 
the reduced equations (I2.5P imply 

G^'^n, = nM^C"^ - ]-n^'VpCP, (2.9) 

where 

riu = d,yt 

is the unit timelike normal to the Cauchy hyperurfaces. Thus if 

G^"",|5o=0, (2.10) 

i.e. if the Hamiltonian and momentum constraints are satisfied by the initial data, and if the reduced equations (j2.4p 
are satisfied then it follows that 

[nM^^^-^n^'VpCP]\s„^0. (2.11) 

It is easy to check that (|2.1ip implies that dtC^\so — provided C^jsj, — 0. 

As a result, the traditional Hamiltonian and momentum constraints on the initial data, along with the reduced 
evolution equations (|2.4p ). imply that the initial conditions (|2.7p and (|2.8p required for preserving the harmonic 
conditions are satisfied. Conversely, if the Hamiltonian and momentum constraints are satisfied initially, then (j2.9p 
ensures that they will be preserved under harmonic evolution. Thus the conditions C^ = can be considered as the 
constraints of the generalized harmonic formulation. 

The formalism also allows constraint adjustments by which (j2.4p is modified by 

Ef"" := C" - V^^C'') + -g'"'S/pCP + Ai^^C = 0, (2.12) 

where the coefficients A^'^ have the dependence A^ (x^ g, dg). Such constraint adjustments have proved to be impor- 
tant in applying constraint damping ^27j in the simulation of black holes ^2 8 - .30] and in suppressing long wavelength 
instabilities in a shifted gauge wave test [3l| (see Sec. IXIli|) . However, they do not change the principal part of the 
reduced equations and have no effect on well-posedness. 

Historically, the first Cauchy codes were based upon the "3+1" or Arnowitt-Deser-Misner (ADM) formulation of 
the Einstein equations [S^l ■ The ADM formulation introduces a Cauchy foliation of space-time by a time coordinate 
t and expresses the 4-dimensional metric as 

ds^ = -a^dt^ + h,j (dx' + P^dt) (dx^ + (3^dt) , (2.13) 

where hij is the induced 3-metric of the t = const foliation, a is the lapse and /?' the shift, with the unit normal to 
the foliation given by n^ = (1, — /3*)/a. 

The field equations are written in first differential form in terms of the extrinsic curvature of the Cauchy foliation 

This can be accomplished in many ways. In one of the earliest schemes proposed for numerical relativity by York [33{, 
the requirement that the 6 spatial components of the Ricci tensor vanish, i.e Rij = 0, yields a set of evolution equations 
for the 3-mctric and extrinsic curvature, 

dtgij - Cpgtj = -2akij (2-14) 

dthj-Cphj = -D,Dja + a{n,j+khj~2k\kij), (2.15) 



where Di is the connection and TZij is the Ricci tensor associated with hij. The Hamiltonian and momentum 
constraints take the form 

2G^"'nf,n^ = 7^ - hjk'^ + fc^ = (2.16) 

C'n^ = Dj {k'^ - h'^k) = 0, (2.17) 

where TZ = h^^TZij and k = h^^kij. 

Codes presently used for the simulation of binary black holes apply the constraints to the initial Cauchy data but do 
not enforce them during the evolution. The choice of evolution equations may be modified by mixing in combinations 
of the constraint equations. In addition, the evolution equations must be supplemented by equations governing the 
lapse and shift. There is a lot of freedom in how all this can be done. The choices affect whether the Cauchy problem 
is well-posed. 

III. COMPLICATIONS OF THE IBVP 



The difficulties underlying the IBVP have recently been discussed in [3J-|37| . There are several chief complications 
which do not arise in the Cauchy problem. 

1 . The first complication stems from a well-known property of the flat-space scalar wave boundary problem 

id^-V^)<^^0, x>0,t>0. (3.1) 

The light rays are the characteristics of the equation. There are two characteristics associated with each direction, 
e.g. the characteristics in the ztx direction. Both of these characteristics cross the initial hypersurface i = but 
only one crosses the boundary at a; = 0. As a result, although the initial Cauchy data consist of the two pieces 
of information $|t=o ^nd dt^\t=o, only half as much boundary data can be freely prescribed at x = 0, e.g the 
Dirichlet data qo — 9t$|x=o, or the Neumann data qjy — dx^\x=o or the Sommerfeld data qs = {dt — dx)'^\x=o- 
Sommerfeld data is based upon the derivative of $ in the characteristic direction determined by the outward 
normal to the boundary. (In the first differential order formalism [dt — dx)'^ is an ingoing variable at the 
boundary.) The choices qo = Qn = Qs do not lead to the same solution. In order to obtain a given physical 
solution, this implies that the boundary data cannot be prescribed before the boundary condition is specified, i.e. 
the boundary data for the solution depends upon the boundary condition, unlike the situation for the Cauchy 
problem. The analogue in the gravitational case is the inability to prescribe both the metric and its normal 
derivative on a timelike boundary, which implies the inability to freely prescribe both the intrinsic 3-mctric of 
the boundary and its extrinsic curvature. This leads to a further complication regarding constraint enforcement 
at the boundary, i.e. the Hamiltonian and momentum constraints (j2.16p and (J2.17I) cannot be enforced directly 
because they couple the metric and its normal derivative. 

2. For computational purposes, a Sommerfeld boundary condition is preferable because it allows numerical noise 
to propagate across the boundary. Thus discretization error can leave the numerical grid, whereas Dirichlet and 
Neumann boundary conditions would reflect the error and trap it in the grid. The Sommerfeld condition on 
a metric component supplies the value of the derivative K°'dag^v hi an outgoing null direction K"' . However, 
the boundary does not pick out a unique outgoing null direction at a given point but, instead, essentially a half 
cone of null directions. This complicates the geometric formulation of a Sommerfeld boundary condition. In 
addition, constraint preservation does not allow free specification of Sommerfeld data for all components of the 
metric, as will be seen later in formulating the Sommerfeld conditions (|10.2ip - (|10.23p . 

3. The correct boundary data for the gravitational field is generally not known except in special cases, e.g. when 
simulating an exact solution. This differs from electromagnetic theory where, say, homogeneous Dirichlet or 
Neumann data for the various components of the electromagnetic field correctly describe the data for reflection 
from a mirror. The tacit assumption in the simulation of an isolated system is that homogeneous Sommerfeld 
data gives rise to minimal back reflection of gravitational waves from the outer boundary. But this is an 
approximation which only becomes exact in the limit of an infinite sized boundary. 

4. Another major complication arises from the gauge freedom. In the evolution of the Cauchy data it is necessary to 
introduce a foliation of the spacetime by Cauchy hypersurfaces St, with unit timelike normal Ua- The evolution 
of the spacetime metric 

gab = -naUb + hab (3.2) 



is carried out along the flow of an evolution vector field t° which is related to the normal by the lapse a and 
shift /3" by 

t" = an" + /3° , /3"na = 0. (3.3) 

The choice of foliation is part of the gauge freedom in the resulting solution but does not enter into the 
specification of the initial data. In the current treatments of the IBVP, the foliation is coupled with the 
formulation of the boundary condition. As a result, some gauge information enters into the boundary condition 
and boundary data. 

5. The partial derivative dagfn, entering into the construction of the boundary condition for the metric has by 
itself no intrinsic geometric interpretation, unless, say, a background connection or a preferred vector field is 
introduced. 

6. In general, the boundary moves with respect to the initial Cauchy hypersurface in the sense that the spacelike 
unit outer normal Na to T is not orthogonal to the timelike unit normal n^ to Sq. The initial velocity of the 
boundary is characterized by the hyperbolic angle 0, where 

TV^ria = sinhe (3.4) 

Specification of 8 on the edge So must be included in the data. 

The coordinate specification of the location of the boundary is pure gauge since it does not determine its 
location geometrically in the sense that a curve is determined geometrically by its its acceleration, given its 
initial position and velocity. Given Bq and 0, the future location of the boundary should be determined in a 
geometrically unique way. In the Friedrich-Nagy system, the motion of the boundary is determined by specifying 
its mean extrinsic curvature. But this is tantamount to a piece of Neumann data. Can this be accomplished via 
a non-refiecting boundary condition of the Sommerfeld type? 

7. In a reduction to first differential order form by introducing a momentum 11, according to the example 

n'^da^ = n, (3.5) 

there is a further difficulty if 7^ at the boundary. The sign of Q determines whether n" points outward or 
inward to T, i.e whether $ is an ingoing or outgoing variable. Thus the sign of G determines whether such an 
advection equation requires a boundary condition. This forces a Dirichlet condition on the normal component 
of the shift in some 3+1 formulations. 

8. There are also compatibility conditions between the initial data and the boundary data at the edge Bf). For 
the example of the scalar wave problem p.ip with a Dirichlet boundary condition at a; = 0, the boundary data 
must satisfy 

a2$|(,^o,.=o) = {dl + dl + af)$|(,=o,.=o), (3.6) 

where the right hand side is determined by the initial data. An infinite sequence of such conditions follow 
from taking time derivatives of the wave equation. They must be satisfied if the solution is required to be 
C°° . In simple problems, this sequence of compatibility conditions can be satisfied by choosing initial data and 
boundary data with support that vanishes in a neighborhood of the edge Bq. But in problems with elliptical 
constraints, such as occur in general relativity, this simple approach is not possible. In numerical relativity, these 
compatibility conditions are usually ignored, with the consequence that some transient junk radiation emanating 
from the edge is generated. In principle, this could be avoided by smoothly gluing the initial data to an exterior 
region with Schwarzschild [38l] or Kerr [39,] data. This gluing construction would avoid mathematical difficulties 
but it is an implicit construction and in practice no numerical algorithm for carrying it out has been proposed. 
In the simulation of binary black holes, this edge effect combines with another source of junk radiation which 
is hidden in the choice of initial data. The tacit assumption is that the spurious radiation from these sources 
is quickly flushed out of the simulation, with no significant effect after a few crossing times. Since this issue is 
difficult to treat or quantify in a useful way, I adopt the expedient assumption that all compatibility conditions 
for a C°° solution have been met. 

In order to resolve most of the above complications it appears that a foliation Bt of the boundary T must be 
specified as part of the boundary data. Such a foliation is a common ingredient of all successful treatments to date. 



8 

The foliation supplies the gauge information which determines a unique outgoing null direction for a Sommerfeld 
condition. In Sec. IIV[ I specify Bt in terms of the choice of an evolution vector field on the boundary. 

Most of the above complications stem from the fact the domain of dependence determined by the boundary alone 
is empty. An initial value problem for a hyperbolic system can be be consistently posed in the absence of a boundary. 
But the opposite is not true for an IBVP. Without an underlying Cauchy problem, an IBVP does not make sense. In 
an IBVP, boundary data cannot determine a unique solution independently of the initial Cauchy data and there is 
no domain in which the solution is independent of the initial data. Thus a well-posed IBVP problem must be based 
upon a well-posed Cauchy problem. 



IV. THE BARE MANIFOLD 



In constructing an evolution code for the gravitational field, the first step is define a spatial grid and a time update 
scheme. This sets up the underlying structure necessary to store the values of the various fields. The analogous object 
at the continuum level corresponds to the bare manifold on which the gravitational field is later painted. This is the 
approach I will adopt here. It provides a useful way to order the introduction of the basic geometric quantities which 
enter the IBVP. 

Setting up the spatial grid corresponds to the analytic specification of spatial coordinates x* on Sq. The time 
update algorithm corresponds to the introduction of an evolution field i" in Ai which is tangent to the boundary 
T. (In more complicated update schemes, which I won't consider, the boundary might move through the grid.) The 
evolution field must have the property that under its flow Sq is mapped into a foliation St of A^ , and its edge Bq is 
mapped into a foliation Bt of T. 

If a time coordinate is initiated at i = on 5o then the flow of i" induces adapted coordinates a;'' — {t, x^) on Jv[ 
by requiring 

Ctt = 1 (4.1) 

Ctx' = 0, (4.2) 

where Ct is the Lie derivative with respect to t°- . Note that t°- and the adapted coordinates x^ are explicitly constructed 
flelds on M. with no metric properties. They uniquely flx the gauge freedom on A^ in precisely the same way that 
the numerical grid and update scheme provide a unique evolution algorithm. If x are the coordinates on the edge 
i3o: then the corresponding adapted coordinates on the boundary are (t, a; ), where Ctx = 0. It will be convenient 
throughout this review to let the manifold with boundary be described by adapted coordinates 

x'' = (t>0,a;> 0,x^). (4.3) 

Under a diffeomorphism ip oi Ai which maps T into itself, i° — ?► ip^f^ where V'*i" can be chosen to be any other 
possible evolution flcld. In particular, if 1p^,t'^ — t"' in A4 and ipx"^ — x"^ on Sq then the diffeomorphism must be the 
identity. Thus, given a coordinate gauge on Sq, the choice of t° determines the remaining diffeomorphism freedom. 
This allows a description of the evolution in a specific choice of adapted coordinates without losing sight of the gauge 
freedom. 

Note that any one-form normal to the boundary is proportional to dat. However, at the bare manifold level the 
unit normal cannot be specified since that involves metric information. The projection tensor 

< = 5t - f'dbt (4.4) 

has the properties 

^tv^'dat = 0, (4.5) 

i.e. it projects a vector field into the tangent space of St, and 

TTtWat'' = 0, (4.6) 

i.e. it projects a 1-form into the space orthogonal to dt- 

These are the main structures that exist in the IBVP a priori to introducing a geometry on A^. There is an 
alternative approach in which geometrical concepts are introduced earlier. In the Cauchy problem, the initial data 



can be specified in geometrical form as tensor fields hat and kab on a "disembodied" 3-manifold Sq (cf. [4^). Only 
after the embedding of Sq in M is this data interpreted as the intrinsic metric hat and extrinsic curvature kab of Sq. 
The mean curvature k = h°'^kab can itself be interpreted as a variable determining the location of iSq. Similarly, in the 
IBVP the mean extrinsic curvature of T can be interpreted as a wave equation determining the geometric location of 
the boundary [161] . However, these interpretations assume knowledge of the spacetime geometry which is only known 
after a solution is found. 

This order in which the basics objects are introduced is akin to the question: Which came first - the geometry or the 
manifold (or some combination)? Here I adopt the manifold approach, which is more akin to the spirit of numerical 
relativity. I assume a priori, for the given choice of evolution field i", that M is the domain of dependence of the 
initial-boundary data, i.e. it is the manifold upon which the data determines a unique evolution. Here "a priori" is 
used in the sense of a spacetime geometry which exists only after the solution of the IBVP is obtained. 



V. INITIAL DATA 

Since Einstein's equations are second differential order in the metric, any evolution scheme must specify 5^1/ and 
dtg^Lv on Sq. The classic result of the Cauchy problem is that a geometrically unique solution of the Cauchy problem 
is determined by initial data consisting of the intrinsic metric hab of Sq and its extrinsic curvature kab^ subject to the 
constraints (P?TCl) - (PTf)) . 

The remaining initial data necessary to specify a unique spacetime metric consist of gauge information, i.e. gauge 
data that affect the resulting solution only by a diffeomorphism. One such quantity is the lapse a, which relates the 
unit future-directed normal to the time foliation according to 

Ua = -adat. (5.1) 

The embedding of Sq in Ai then gives rise to the spacetime metric 

gab = -naUb + hab (5.2) 

and the interpretation of kab as the extrinsic curvature through the identification 

kab = K^c-rib, (5.3) 

where Vc is the covariant derivative associated with gab- 

The choice of evolution field t° supplies the remaining gauge data. It is transverse but not in general normal to the 
Cauchy hypersurface so that it determines a shift /3a according to 

/3a - habt''. (5.4) 

This relationship supplies the metric information 

gabt'' = aUa + Pa 

relating i° to the unit normal Ua- 

In the adapted coordinates, the metric has components 

gtt = -a' + h^jP'P' (5.5) 

gu = /3^ = h^jPi (5.6) 

gtj = h^j. (5.7) 

The inverse metric is given by g"'' = —n°'-n!' -t- /i"'', where 

h^^Ub = 0, h'^'hbc = 51+ n^Ub. (5.8) 

In the adapted coordinates, 

5" = -a-' (5.9) 
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g*^ = a-'^P' (5.10) 

g^' = h'^ , h^'^kkj = 5). (5.11) 

The implementation of the initial data into an evolution scheme depends upon the details by which Einstein's 
equations are converted into a set of PDEs governing g^^, in the adapted coordinates. All such schemes require 
specification of the initial values of the lapse and shift, in addition to hij and fc^. Thus it can be assumed that 
gab is specified on 5o. By Lie transport along the streamlines of i", this then allows the construction of a preferred 
stationary background metric g^f^ on M picked out by the initial data. Given the choice of evolution field t°- and the 
initial Cauchy data, this background metric is uniquely and geometrically determined by 

^tgab = 0> 9ab\So = 9ab\so- (5-12) 

In the adapted coordinates g^i^it^x'') — 17^1,(0,0;*). 



VI. HYPERBOLIC INITIAL-BOUNDARY VALUE PROBLEMS 



There is an extensive mathematical literature on the IBVP for hyperbolic systems. The major progress traces 
back to the formulation of maximally dissipative boundary conditions for linear symmetric hyperbolic systems due to 
Fricdrichs [41| and Lax and Phillips 42] . There has been recent progress in obtaining results for quasilinear systems 
where the boundary contains characteristics, as arises in some formulations of Einstein's equations. Unfortunately, 
much of this material is heavy on the mathematical side and not easy reading for relativists coming from astrophysical 
or numerical backgrounds. In the absence of the complications of shocks introduced by hydrodynamic sources, 
relativists are content to deal with smooth, i.e. C°°, solutions and forgo the Sobolev theory which enters a complete 
discussion of the quasilinear IBVP. For relativists, the most readable source on the theory of hyperbolic boundary 
problems is the textbook by Kreiss and Lorenz [43], which boasts: "In parts, our approach to the subject is low- 
tech.... Functional analytical prerequisites are kept to a minimum. What we need in terms of Sobolev inequalities is 
developed in an appendix." Taylor's 44, 45] treatises on partial differential equations contain a classic treatment of 
pseudo-differential theory but are less readable for relativists. Fortunately, much of the critical formalism pertinent 
to Einstein's equations appears as background material in papers on the gravitational IBVP. The material I present 
here is heavily based upon those sources, namely |15l - [l9l . [37ll43j . 

There are two distinct formulations of the IBVP, depending upon whether you consider Einstein's equations as a 
natural second differential order system of wave equations or whether you reduce it to a first order system. While 
the second order approach is the most economical, it is not applicable to all formulations of Einstein's equations, 
particularly those whose gauge conditions do not have the semblance of wave equations. The first order theory has been 
extensively developed because of its historic importance to the symmetric hyperbolic formulation of hydrodynamics. 
The IBVP for second order systems has received less attention and some new techniques have originated in the 
consideration of the Einstein problem. 

There are also two distinct approaches to studying well-posedness - one based upon energy estimates and the other 
based upon pseudo-differential theory where Fourier-Laplace expansions are used to reduce the differential operators 
to algebraic operators. The pseudo-differential theory can be applied equally well to first or second order systems. In 
the following, I give a brief account of the underlying ideas in terms of some simple model problems. This will provide 
a background for discussing the difficulties that arise when considering constraint preservation in the gravitational 
IBVP. 

The subclasses of hyperbolic systems consist of weak hyperbolicty, strong hyperbolicity, symmetric hyperbolicity 
and strict hyperbolicity. These subclasses are determined by the principal part of the system when written as first 
differential order PDEs. Weakly hyperbolic systems do not have a well-posed Cauchy problem, which turned out to be 
responsible for the instabilities encountered in early attempts at numerical relativity using naive ADM formulations 
with a prescribed lapse and shift. Strong hyperbolicity is sufficient to guarantee a well-posed Cauchy problem but 
not a well-posed IBVP. It is possible to base a well-posed IBVP on either symmetric hyperbolic or strictly hyperbolic 
systems. However, strictly hyperbolic systems arise very rarely and, to my knowledge, not at all in numerical relativity. 
So I will limit my discussion to the symmetric hyperbolic case. 

I begin with the IBVP for a second order scalar wave equation, where the underlying techniques are transparent 
rather than hidden in the machinery of symmetric hyperbolic theory. The generalization to systems of quasilinear wave 
equations, described in Sec. IVI A31 can also be treated by the same techniques as for symmetric hyperbolic systems. 
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However, for application to the harmonic Einstein problem, the scalar treatment suffices since the principal part of 
the system consists of a common wave operator acting on the individual components of the metric. The mathematical 
analysis which is necessary for a treatment of the quasilinear IBVP in full rigor is beyond my competence and 
presumably outside the interest of someone from a more physical or computational background. I simply state the 
main results and give references when such mathematical theory must be evoked. 



A. Second order wave equations 



The ideas underlying the well-posedness of the IBVP are well illustrated by the case of the quasilinear wave equation. 
I give some examples which are are relevant to the harmonic formulation of Einstein's equations and which illustrate 
the techniques behind both the energy approach and the pseudo-differential approach. 



1. The energy method for second order wave equations 
First consider first the linear wave equation for a scalar field, 

'^tt=^xx+'^yy + ^zz+F (6.1) 

on the half-space 

x > 0, — cx) < y < oo, — oo < z < oo, (6.2) 

with boundary condition at a; = 

$t - a$, - /32$j, - /33$. = Q , a>0, (3l+(3l< 1, (6.3) 

with boundary data q, initial data of compact support 

* = /i, $t=/2, i = (6.4) 

and forcing term F{t, x, y, z). The subscripts {t, x, y, z) denote partial derivatives, e.g 

at 
All coefficients and data are assumed to be real and a > 0, P2, Pz are constants. The notation 

($,*), |i$||2 = ($,$); ($,Vl/)s, ||$||| = ($,<i>)B, 

is used to denote the L2-scalar product and norm over the half-space and boundary space, respectively. 

In order to adapt the standard definition of energy estimates to second order systems, the notation # = 
($, $t, $2:,<I>j^, $2) is used to represent the solution and its derivatives; and similarly f = {fi, f2, fix, fiy, fiz) for 
the initial data. 

For the scalar IBVP (|6.ip - (|6.4p . strong well-posedness requires the existence of a unique solution satisfying the a 
priori estimate 

\mt)r + £ ||*(r)|||dT < Kr (\\ir + f^ WnrWdr + f^ \\q{T)\\ldT^ , (6.5) 

in any time interval < i < T, where the constant Kt is independent of F , f and q. 

It is important to note that (|6.5p estimates the derivatives of 4>, both in the interior and on the boundary, in terms 
of the data and the forcing. This is referred to as "gaining a derivative" . This property is crucial in extending the local 
IBVP to global situations, e.g. where the boundary is a sphere or where there is an interior and exterior boundary 
as in a strip problem. Otherwise, refiection from the boundary could lead to the "loss of a derivative" , which would 
lead to unstable behavior under multiple reflections. 
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The usual procedure is to derive an energy estimate by integration by parts, using for example ($^^,$2) = 
— {^y,^yz) = 0. Consider first the estimates of the derivatives of $ in the homogeneous case F = q = 0. Using 
the standard energy for a scalar field, integration by parts gives 

dtmtf + \\<i>.f + \\%f + \\<^yf) = -2($t, $,)s. 

li (32 = Pa — and a > in the boundary condition (|6.3p then ($t, ^x)b ^ 0, i.e. the boundary condition is dissipative, 
and there is an energy estimate. Otherwise there is no obvious way to estimate the boundary flux. Instead, it is 
possible to use a non-standard energy E for the scalar wave equation (|6.1|) which does provide the key estimate if 
/3| + I3i > 0. 

The first step is to show that 

E := ||$J2 + ||$^||2 ^ ||^^||2 ^ ||^^||2 _ 2($t,/32$j, + /33$,) (6.6) 

is a norm for the derivatives {^t,^x,^y,^z)- Since Z?! + /^l < 1: this follows, after a rotation, from the inequality 
$2 + ^2 _ 2^$t* > for ^2 < 1. 

This leads to 

Lemma 1: The solution of (|6.ip - (|6.4p satisfies the energy estimate 

d,E + a\\^x\\l<E+\\Fr + -U\l- (6.7) 

a 

Proof. Integration by parts gives 

dtW^tf = 2($t, $tO = -dt{\\<^xf + \\<^yf + W^zf) + 2i<Pt,F) - 2($t,$,)B (6.8) 

and 

2dti^t,l32'^y + /33*.) = 2($tt, 132% + /33$.) = -2($„ /32$„ + /33$.)b + 2(F, (32% + /?3*.). (6.9) 

Since (|6.3p implies 

2(«>t, $:.)b = 2a|l$^|i| + 2($^, /32$y + /Ss^^)^ + 2($:„ g)ij, 

subtraction of (|6.9p from (16. 8p leads to 

a*^ = 2($t-/32$j,-/33$.,F)-2a||$,|||-2($„<z)s 

< ||$t - P2% - /33$.f + \\Ff - a\\%\\l + i||<z|||. 



The identity 



\% - (32% - /33$.f ^E- \\%r - \\%r - W'i'.r + ii/?2$. 



^y-PS-^zW ^ ^ - W'i'xW - W-i-yW -W-^zW + || P2'i'y "t" ^3^^ II 

then implies (j6.7p and proves the lemma. 
By integration, the lemma estimates 

^(T) and / \\<i>a:\\%dt in terms of £'(0), / Hi^H^di and / \\q\\%dt. 
Jo Jo Jo 

Strong well-posedness (|6.5p also requires estimates of the boundary norms H^fUs, ||$j^||b and ||$i||_B. First, a 
calculation similar to that above gives 

dti^x,%) = {%U%) + {'^cc,%t) 

= ^l\\%\\l -l\\%\\l + l\\%\\l + ^ll^.lll + i%,F)- (6.10) 
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Estimates of ||<&t||B in terms of |1$j,|1_b, H^^Hb, H^kHb and H^lls can be obtained from the boundary conditions (16.31) . 
which give, for any 6 with < S < 1, 

ll'Stlll = \\(32'^y + l33^z+a<i>^ + q\\l 

< \\l32^y + P3^z\\l + 2\\(32%+l33<^,\\B\\a<^x + q\\B + \\a^x + q\\l (6.11) 

< (1 + S)\\P2'^y + Pa^zWl + (1 + ^)ll«*:. + gill 

< (1 + 6) {Pi + /3|) (||$,||| + ||$,|||) + (1 + i)||a$, + q\\l. 

Next, since /3| + /3| < 1, 5 can be chosen such that (1 + (5)(/3| + /3|) < (1 - S). Therefore, by (plll| . 

^mvWl + W^zWl) < (1 + ^)ll"*:. + q\\l + ll*.|ll + 29t(<i>.,$t) - 2(<i>„F). 

Since ($a;, $() can be estimated by E, there follows 
Lemma 2: 



<const{E{0)+ \\Ffdt+ \\q\\%dt). 

Jo Jo 



An estimate for $ itself can easily be obtained by the change of variable $ — >■ 6***$, as described in Sec. IVI A2l or 
in Appendix 1 of [1^. Together with the results of Lemma 1 and Lemma 2, this establishes 



Theorem 1: The IBVP (j6.ip - (|6.4p is strongly well-posed in the sense of (|6.5p . 

The result can also be generalized to half-plane problems for wave equations of the general constant coefficient form 

$ti = 2b'^,t + /I'-'^jj, xi > 0, -cx) < {y, z) < cx), a;* == (a;, y, z) (6.12) 

where W^ is is a metric of (+ + +) signature. By coordinate transformation, (|6.12p can be transformed into (|6.ip and 
the appropriate boundary conditions formulated. See [18| for details. 

This example illustrates from the PDE perspective the constructions necessary to establish strong well-posedness. 
For the purpose of establishing the strong well-posedness of the IBVP for the wave equation on a general curved space 
background, it is also instructive to take advantage of the geometric nature of the problem. 

In terms of standard relativistic notation, consider the wave equation 

g-^^VaVb* - F (6.13) 

for a massless scalar field propagating on a Lorentzian spacetime M. foliated by compact, 3-dimensional time-slices 
iSt, with boundary T foliated by Bt- Here V denotes the covariant derivative associated with the spacetime metric 
g°-^ . For notational simplicity, let $a = Va$. 

The IBVP consists in finding solutions of ()6.13p subject to the initial Cauchy data 

*l5o=/i. ri'^^,\^^^^h (6.14) 

and the boundary condition 

[(T'' + aiV^)$6]^ = g, (6.15) 

with data q on T ■ Here n is the future-directed unit normal to the time-slices St and and N is the outward unit 
normal to T; T^ is an arbitrary future-directed timelike vector field which is tangent to T; and a > 0. The motion of 
the boundary is described geometrically by the hyperbolic angle tanhG = N^ni,. Without loss of generality, assume 
the normalization g^cT^T'^ = — 1. A Sommerfeld boundary condition then corresponds to the choice a — 1 for which 
T'' + N^ points in an outgoing null direction. 

In order to establish estimates, consider the energy momentum tensor of the scalar field 

8^ = $.$" - i(5?$=$c- 
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The essential idea is the use of an energy associated with a timehke vector m" = T° + SN"", where < 5 < 1, so that 
w° points outward from the boundary. The corresponding energy E(t) and the energy flux P{t) through Bt are 



's, 
which is a covariant version of the non-standard energy (j6.6p , and 

,bc\a 1 



E{t) = / v'Qlna, (6.16) 

J St 



P{t) = j u'QlNa. (6.17) 

Bt 



SO that 



It follows from the timelike property of u"" that E{t) is a norm for $a(0- 
Energy conservation for the scalar field, i.e. integration by parts, gives 

dtE = P- ({QabV'y!' + u'^^aF), 

St 

dtE <P + const{E + j F^). (6.18) 

St 

The required estimates arise from an identity satisfied by the flux density 

u''Q^,Na = -- {{N'^'Paf + {T^-^af + Q°'$a$&) + N'^'^'aT^^b + <5(7V"$a)' + S{T''<^af (6.19) 

where Qbc = 9bc + TbTc — N^Nc is the positive deflnite 2-metric in the tangent space of the boundary orthogonal to 
T'^. By using the boundary condition to eliminate r°$a in the last group of terms, there follows 



u 



^e'^Na = --((iV"$„)2 + (T''$,)2 + Q''&$,$,) 

+ (-a + (5(1 + a^)) (iV°$a)^ + (1 - 2aS)N''<i>aq + Sq^ 

= -^((iV'^$,)2+(T'^$,)2+Q''b$„$,) 



2 

2(a-5(l-f-a2))J ' V" ' 4(a-(5(l + a2)) 



(»-.(l + a-))(N-.„- ,''-f,°l'' ) +f^+ , ''-f,°i'.J ,' 



The choice 



0<S< 



l + a2 
(which also guarantees that (5 < 1 so that u° is timelike), gives the inequality 



u 



O'iNa < ~- ((A^"$a)' + (T'^^a)' + 0"'$a$;,) + const q\ (6.20) 



2 



It now follows from ((08t and (lOO)) that 

dtE + I ^('(iV'^$,)2 + (r'^$^)2+Qa6^_^$^ 

< const(E+ I F^+ I qA. (6.21) 

This is the required estimate of the gradient $a on the boundary (as well as usual estimate of the energy E) to prove 
that the problem is strongly well-posed. As in the previous example, an estimate of <& itself follows from the change 
of variable $ — > e^*$, which introduces a mass term in (|6.13p . 
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2. The quasilinear case 

The estimate (I6.21[) is sufficient to establish the criterion (I6.5P for strong well-posedness of the IBVP for the hnear 
wave equation with constant metric coefficients. In order to extend the result to the quasilinear case on a curved 
space background, where the metric depends upon $ and $a, it is necessary to show that the corresponding estimates 
hold for arbitrarily high derivatives of $. In the process, this also requires stability of the system under the addition 
of lower differential order terms, which arise under the differentiation of the wave equation. These requirements are 
sometimes neglected or misunderstood in the relativity literature. 

More generally, local existence theorems for variable coefficient or quasilinear equations follow by iteration of 
solutions of the linearized equations with frozen coefficients. The energy estimates for the frozen coefficient problem 
establish the existence of a unique solution which depends continuously on the data. The extension of this result 
to the quasilinear case first requires that the problem with variable coefficients be strongly well-posed. For this, it 
already necessary to obtain estimates for arbitrarily high derivatives of the solution to the linearized problem. 

For the purpose of illustrating the procedure, it suffices to consider the IBVP for the 2(spatial)-dimensional wave 
equation with variable coefficients 

$tt = P$ + i?$ + F, a;>0, -oo<j/<oo, (6.22) 

with smooth initial data 

$ = /i, $t = /2, i = 0, (6.23) 

and boundary condition 

a(t, y)$t = <i>a; - ^<I> + r(t, j/)$ + 5(i,y) , a{t,y) > const > , ^ ^ const > (6-24) 

with smooth, compatible boundary data q. Here 

P(f> == {a<^x)x + (6$y)y - 2/i<i>t - /i^<I> , a > flo = const > 0, b > bo = const > , 

is an elliptic operator which has been modified by terms which arise in (|6.22[) by the transformation <I> — ;> e''*<I', where 
/i introduces a mass term; and 

are terms of lower (zeroth and first) differential order. The coefficients a, b and Ci are smooth functions of {t,x,y). 
Consider the norm for the Cauchy data 

E = ||<i>,||2 + ($„ aa>,) + {%,b<^y) + /i2||$||2. 

Integration by parts leads to 

dtE = ~M\^tf + 2i<Pt,F) + 2($t, m) - 2($t, a$^)B + (at$^, $^) + {bt%,%) 

< const{\\F\\^ + E)-2{<^t,a<^,)B. (6.25) 

The boundary condition gives 

-{<^t,a<^x)B ^ -i<^t,aa^t)B~ K'^t,a<^)B + {'^t,ar<^ + aq)B 

^ -($t,aa$t)B-M(*t,ao*)s-M(*t,(a-ao)$)B + (*t,ar$ + ag)i3 (6.26) 

< -^fiaodt\ml~^{^uaa<Pt)B+const{ml + \\q\\l). (6.27) 



Therefore, from (lOS]) . 



dt{E + fiaoWml) + {^t,aa<^t)B (6-28) 

< const {E + ||<D||| + ||F||2 + ||g|||) , (6.29) 



This establishes that the energy estimate is stable against lower order perturbations. Now it is possible to estimate 
the derivatives. For the derivatives Y — ^y and T — ^t tangential to the boundary, differentiation of the wave 
equation gives 

Ytt = PY + RY + Ry'i + iay'i>,), + ibyY)y + Fy (6.30) 

Ttt = PT + RT + Rt^ + (at'S?,), + {btY)y + Ft. (6.31) 
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Here Ry^ and Rt^ are linear combinations of first derivatives of $ which have already been estimated and can be 
considered part of the forcing term F. Also, the wave equation (16.221) implies 

a^xx = Tt — bYy + terms that have already been estimated, 

so that ^xx is also lower order with respect to (j6.30p and (I6.31|) . Thus, except for lower order terms, Y and T solve the 
same wave equations (j6.30p and (j6.3ip as $, with the same boundary conditions up to lower order terms. Therefore 
all second derivatives ^ay and ^at can be estimated as well as $xx- Repetition of this process gives estimates for any 
number of derivatives. 

In order now to show the existence of a solution to the variable coefficient problem, one approach, which is par- 
ticularly familiar to numerical relativists, is to approximate the PDE by a stable finite difference approximation. 
This approach is detailed in [43| where summation by parts (SBP) is applied to the finite difference problem in the 
analogous way that integration by parts is used above in the analytic problem. The SBP approach shows that the 
corresponding estimates hold independently of gridsize. Existence of a solution of the analytic problem then follows 
from the limit of vanishing gridsize. 

It also follows from the estimates for arbitrary derivatives of the variable coefficient problem that Sobolev's theorems 
can be used to establish similar, although local in time, estimates for quasilinear systems. Then the same iterative 
methods used for first order symmetric hyperbolic systems can be used to show that well-posedness extends locally in 
time to the quasilinear case. In this way it was shown in [1^ that the general quasilinear wave problem (|6.13l) - (|6.15p . 
where the spacetime metric now depends upon $ and $0, is strongly well-posed. 



3. Systems of wave equations 

The strong well-posedness of the IBVP for the quasilinear scalar wave (|6.13p ~ (|6.15p can be generalized to a system 
of coupled wave equations 

(7''^($)VaVf,$^ = i^^($,V$), A = 1,2,...N (6.32) 

with smooth initial data 

^X^^ff, n'V.^^l^^^f^, (6.33) 



and the boundary condition 

(T^ + aN'') Vb$^|^ = c"^B Va$^|^ + d^B $^|^ + q^, (6.34) 

where a — a(a;, $) > 0, g'^ = q^{x), c'^'^b — c'^^b{x,^) and d^B — d'^Bix,^) are smooth functions of their 
arguments. All data are compatible. As before, each time-slice St is spacelike, with future-directed unit normal 
n°'{x, $), the boundary T is timelike with outward unit normal N°'{x, $) and T" ~ T"-{x, $) is an arbitrary future- 
directed timelike vector field tangent to T. 

In [1^, it was shown that the IBVP (|6.32p ~ (|6.34[) is strongly well posed given certain restrictions on c'^^b, i-e. 
there exists a solution locally in time which satisfies (|6.5p in terms of the corresponding L2 norms for ^^ and its 
gradient. One important situation in which the restrictions on c" b are satisfied is when it can be transformed into 
the upper diagonal form 

c'^^s^O, B<A. 

This has important applications to the constrained systems of wave equations obtained by formulating Maxwell's 
equations in the Lorentz gauge or Einstein's equations in the harmonic gauge [19|, as discussed in Sec. 1X1 



4. Pseudo-differential theory 



The previous scalar wave problems required a non-standard energy to obtain the necessary estimates. In more 
general problems, an effective choice of energy might not be obvious or even exist. Pseudo-differential theory provides 
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an alternative treatment of such cases. The approach is based upon a Fourier transform for the spatial dependence 
and a Laplace transform in time. It can be applied equally well to first or second order systems. Here I illustrate how 
it applies to the second order wave equation. The more general theory for a system of equations is usually presented 
in first order form and is reviewed in Sec. IVIB II 

As an illustrative example, consider the 2(spatial)-dimensional version of the IBVP for the scalar wave considered 
in Sec. IVIAll 

'^tt = ^xx + ^yy + F, x>0, ~oo<y< cx), (6.35) 

with boundary condition at a; = 

$t - a$, - /3$.y = g, a > 0, |/3| < 1, (6.36) 

with compact boundary data q and initial data 

^0,x,y) = h{x,y), M0,x,y)^f2{x,y). (6.37) 

The following simple observation reveals the underlying idea. The system cannot be well-posed if the homogeneous 
version of (|6.35p - (|6.37l) with F = q = admits arbitrarily fast growing solutions. The homogeneous system has 
solutions of the form 

$(i, X, t) = e^*+^"V(^), IV'U < oo, (6-38) 

where 

'Pxx - (s^+w^)^ = 0, (6.39) 

s^(0) = a(^:r(0)+i/3a;^(0). (6.40) 

Here ^p{x) is a smooth bounded function, so that its maximum norm |(p|oo is finite, w is a real constant and s = rj + i^ 
a complex constant. This poses an eigenvalue problem. If there are solutions with r; = Jfts > then 

$^ = e'^^'+^^^y ip{vx) 

is also a solution for any i^ > 0, so that there are solutions which grow arbitrarily fast exponentially. Therefore, 
a necessary condition for a well-posed problem is that solutions with 5R s > must be ruled out by the boundary 
condition. 

The general solution of the ordinary differential equation (|6.39p is 

^ = aie''+^+a2e''-^, (6.41) 

where 



K±=±^Js^+^ (6.42) 

solve the characteristic equation 

K^ - (S^+W^) =0. 

Here 3?^+ > and 3fiK_ < for 3fis > 0. By assumption ^p is bounded which requires that ci = 0, i.e. 

ip^a^e''-''. (6.43) 

Introducing (|6.43p into the boundary condition (|6.40p gives 

(s - aK_ - il3uj)a2 = 0. (6.44) 

Since 3fiK„ < and by assumption a > 0, there are no solutions with IKs > 0. Thus this necessary condition for a 

well-posed problem is satisfied. 

In order to proceed further it is technically convenient to assume that the initial data (j6.37p vanishes. This may 
always be achieved by the transformation 

$^$-e-Vi-te-*(/2 + /i), (6.45) 
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so that the initial data gets swept into the forcing term F. Then (|6.35p - (|6.37l) can be solved by Fourier transform 
with respect to y and Laplace transform with respect to t, i.e. in terms of 

-I /"OO /"OO 

^{s,x,uj)^^= dye-"^y dte-''^t,x,y), 3fis>0, (6.46) 

V27r J-oo Jo 

where lo is real and s is complex. The inhomogeneous versions of (j6.39p and (j6.40p imply that the coefficients satisfy 

(s - i;3w)$(s,0,w)-a$^(s,0,w) =q(s,w). (6.47) 

Since it has already been shown that the homogeneous system (|6.39p - (|6.40p has no eigenvalues for 3ff s > and 
5ft K < 0, it follows that (j6.47p has a unique solution $. Inversion of the Fourier-Laplace transform then gives a unique 
solution for 4>. 

The well-posedness of a variable coefficient or quasilinear wave problem also requires estimates of the higher deriva- 
tives of <I>. The system of equations for the derivatives are obtained by differentiating the wave equation and the 
boundary condition. In that process, any variable coefficient terms in the boundary condition lead to inhomogeneous 
boundary data for the derivatives. It is possible to transform the boundary data g to by a transformation anal- 
ogous to (J6.45I) . which sweeps q and its derivatives into the forcing term F. If this inhomogeneous boundary data 
is continually subtracted out of the boundary condition, inhomogeneous terms of higher differential order appear in 
the forcing term. As a consequence, the resulting estimates would bound lower derivatives of the solution in terms of 
higher derivatives of the data, a process referred to as "losing" derivatives. 

Instead, a different approach is necessary to establish well-posedness. It is simple to calculate the solution for 
F = 0. Corresponding to (|6.43p and (|6.44p . there follows 



where 



$(s,a;,w) =e''-^$(s,0,cj) (6.48) 

(s — aK_ — i/3a;)$(s, 0, a;) — q{s, uj). 
It is now possible to establish |17| 
Boundary stability: The solution (|6.48p satisfies the estimates 

|$,(s,0,w,)| < K\q{s,Lj)\, 

^/W+^-ms,0,Lo)\ < K\q{s,Lu)l (6.49) 

where the constant K is independent of s and uj. Similar estimates hold for all the derivatives. 

The estimates (I6.49P follow from purely algebraic consequences of the eigenvalue relations (|6.42p and (|6.44l) . The 
essential steps are to show 

1. There is a constant Si > such that|5RK_| > Sidis 



2. For all w and s with 5Rs > 0, there is a constant ^2 > such that \s — an^ — i(3uj\ > (52^|sp + |wp. 
See [l3l for the details. 

Boundary stability allows the application of the theory of pseudo-differential operators to show that the IBVP is 
strongly well-posed in the generalized sense, 



/*||*(T)fdr+ f\mT)\\ldT<KT(f\\F{T)fdT+ f\\qiT)\\ldr), 0<t<T, 
Jo Jo \Jo Jo / 



(6.50) 



where * = ($, $a)- 

The theory is discussed in Sec. IVIB II in the standard context of first order systems. But the first order theory is 
flexible enough to apply to second order systems. In particular, in Sec lVlBTj it is applied to show that the IBVP for 
the quasilinear version of the second order wave equation with boundary condition (j6.35p - (|6.36l) is well-posed in the 
generalized sense. 

Strong well-posedness in the generalized sense is similar to strong well-posedness ()6.5p except now the initial data 
has been swept into the forcing term and the estimate for $ in the interior involves a time integral. In both cases, 
the gradients at the boundary are estimated by the boundary data and the forcing, i.e. a derivative is gained at the 
boundary. This ensures that the well-posedness of the local halfplane problem can be extended globally to include 
boundaries that lead to multiple reflections. 
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5. Generalized eigenvalues 

Strong well-posedness in the generalized sense not only rules out eigenvalues of (I6.40p with rj — ^{s) > but also 
generalized eigenvalues for which r; = 0. This is implicit in the estimates for boundary stability (I6.49P in which the 
constant K is independent of s. However, generalized eigenvalues can exist in well-behaved physical systems. A prime 
example is a surface wave which travels tangential to the boundary with periodic time dependence. See [l5| for the 
treatment of such an example from Maxwell theory. 

Generalized eigenvalues are ruled out by the boundary conditions required for strong well-posedness in the gener- 
alized sense and historically have been treated on an individual basis. However, a new approach to this problem has 
recently been formulated by H-0. Kreiss [47|. This approach splits the problem into two subproblems: 

1. One in which the forcing vanishes, _F = 0. 

2. One in which the boundary data is homogeneous, q = 0. 

A second order wave problem is called well-posed in the generalized sense if these subproblems satisfy the corre- 
sponding estimates 

/ mr)\\ldT <Kt [ \\q{T)\\ldT, 0<t<T, (6.51) 

Jo Jo 

( \\^{T)fdT<KT I \\F{T)fdT, 0<t<T. (6.52) 

Jo Ja 

Here it is only required that $, and not its gradient $, be estimated by the boundary data q. Thus the solution 
no longer gains a derivative at the boundary. However, no global problems arise from multiple reflections because the 
estimate (|6.52p implies the gain of one derivative in the interior. 

As examples, consider the scalar wave problem (|6.35p with the two choices of boundary conditions at x = 0, 

(A) $^ - ij3^y = g, 

(B) $, _ p^y ^ g, 

where /3 is real, with |/3| < 1. In case (A), $ is complex. Introducing these boundary conditions into the homogeneous 
system ((O8| - ((O0l) gives 

(A) K = ujP, 

(B) K = -iujp, 

with s^ — k"^ — uj"^. 

In neither case is there a solution with SRs > but both cases possess generalized eigenvalues, 

(A) s^ = -uj^{l~l3^), SRs^O, 

(B) s^ = -uj^{l+p^), 3?s = 0. 

The corresponding eigenfunctions are the surface waves 



(A) e"^i±^A^t+y)~\-'P\^^ (6.53) 

and the oscillatory waves 



(B) g^^{±V^Tp^^+P-+y) ^ (6.54) 

which give rise to glancing waves for (3 — 0. By investigating the inhomogeneous problem, it can be shown that 
both choices of boundary condition give rise to an IBVP which satisfies (|6.5ip and (|6.52l) and is well-posed in the 
generalized sense [47| . 
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B. First order symmetric hyperbolic systems 

Most of the work on the IBVP for hyperbohc systems has been directed toward fluid dynamics, where a first order 
formulation is natural. Here I describe the essentials of the two main approaches, pseudo-differential theory and the 
theory of symmetric hyperbolic systems. 



1. Pseudo-differential theory 

In order to summarize the pseudo-differential theory for a first order symmetric hyperbolic system consider first 
the constant coefficient system 

m 

ut = V{d^)u + F, V{d^) = V'd^^ =Ad^,+Y, Bjd^, (6.55) 

i=2 



on the half-space 



with initial data 



i > 0, xi > 0, — oo < Xj < GO , j — 2, . . . ,m, 



u{0,x)^f{x). (6.56) 

Here u{t,x) = (u^^'{t,x), . . . ,u^^'{t,x)) is a vector valued function of the real variables {t,x) = {t,xi,. . . ,Xm) and 
A, Bj are constant N x N matrices. In applications to spacetime, m = 3 but the number of spatial dimensions 
does not complicate the theory. The notation {u,v) and |up = {u,u) denotes the inner product and norm in the 
iV-dimensional linear space. All data are smooth, compatible and have compact support. 
The symbol representing the principal part of the system, 

rn 

V{iuj) = iAuji + iB{uj_), B{uj_) = ^ BjUjj \uj\ = 1, (6.57) 

is obtained by replacing dx by its Fourier representation iui = (iwi, iw-), iv^ — {uj2, ■ ■ ■ , ^m)- Symmetric hyperbolicity 
requires that A and B be self-adjoint matrices so that the eigenvectors oiV form a complete set with purely imaginary 
eigenvalues for all real lo. More precisely, there exists a symmetric, positive definite symmetrizer H such that HA 
and HBj are self-adjoint. See |48| for more general applications. 

Here it is also assumed that the boundary matrix A is nonsingular so that it can be transformed into the form 

^=['f Ah)^ (6-58) 

where A^, A^^ are real positive definite diagonal matrices acting on the P dimensional subspace u^ and the {N — P) 
dimensional subspace u , respectively. The theory also applies to the singular case where the boundary is uniformly 
characteristic, i.e. the kernel of A has constant dimension ^3|- See Sec. IVI 1^2] for a treatment of the singular case by 
the energy method. 

The IBVP requires P boundary conditions at xi ~ 0, corresponding to the P ingoing modes in the plane wave 
decomposition carried out below in conjunction with (|6.6ip and (|6.64p . They are prescribed in the form 

u{t,0,x^) = Su (t, 0,x_) + g(i,x_), x- — {x2, ■ ■ ■ ,Xm). (6.59) 

The main ingredient of a definition of well-posedness is the estimate of the solution in terms of the data. See 
Sec. 7.3 of [431. ^Y ^ transformation analogous to (|6.45p . the IBVP (I6.55|) . (|6.56p . (j6.59p reduces to a problem with 
homogeneous initial data / = 0. The required estimate is the first order version of the estimate (j6.50p for the second 
order wave equation. The first order problem is called strongly well-posed in the generalized sense if there is a unique 
solution u such that 

f MrWdr + f \\u{T)\\ldT <kJ f WFirWdT + f |k(r)|||dr) , < t < T, (6.60) 

Jo Jo \Jo Jo J 
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where the constant Kt does not depend on F or q. Here |lw|| and |1w|1_b are the L2 norms of |u| over the half-space 
and the boundary, respectively. 

As for the scalar wave problem considered in Sec. IVI A 4( the IBVP is not well-posed if the homogeneous system 
F = q = Q admits wave solutions 

rri 

M(t,a;i,x_) = e''*+*<"'^>-(/3(xi), (a;,x)_ = ^WjX^, 3fis>0, (6.61) 

with |(/3|oo < 00 • The existence of such homogeneous solutions would imply the existence of solutions which grow 
arbitrarily fast exponentially. 

In order to decide whether such homogeneous waves exist, introduce (|6.6ip into (|6.55p and (|6.59l) to obtain 

Sip = Aipxx+iB{u)-)ip, xi > 0, 
/(O) = ^(^"(0), |(^U < 00. (6.62) 

This is an eigenvalue problem for a system of ordinary differential equations which can be solved in the usual way. 
Let K denote the solutions of the characteristic equation 

Det|AK-(s/-iB(a;_))| =0, (6.63) 

obtained by setting (^(si) = e'^'-^^if{Q). 
It can be shown that 

1. There are exactly r eigenvalues with SRk < and n — r eigenvalues with SRk > 0. 

2. There is a constant ^ > such that, for all s with real 5Rs > and all w_, 

\'Rk\ >5^s. 
In particular, for 5R s > 0, there are no k with ^tn — Q. 
See [4g| for the proof. 

If all eigenvalues Kj are distinct, the general solution of (|6.62p has the form 

where hj are the corresponding eigenvectors. (If the eigenvalues are degenerate, the usual modifications apply.) For 
bounded solutions, all Uj in the second term are zero. Introducing if into the boundary conditions (|6.62p at xi = 
gives a linear system of r equations for the r unknowns (cti, . . . ,ar) = a; of the form 

C(s,w_)o:==0. (6.65) 

Therefore the problem is not well-posed if for some w_ there is an eigenvalue sq with 5R so > 0, i.e. Det C(so, cj_) = 0. 
In that case, the linear system (|6.65p . and therefore also (|6.62p . has a nontrivial solution. Thus the determinant 
condition 

DetC(so,w_) 7^0, 5Rso>0, (6.66) 

is necessary for a well-posed problem. 

Thus in order to consider a well-posed problem assume that Det C 7^ for JR s > 0. Then the inhomogeneous IBVP 
(j6.55p . (J6.56I) . (J6.59I) can be solved by Laplace transform in time and Fourier transform in the tangential variables. 
Again set u(0, x) = f{x) = 0. Then 

su = Aux+iB{uj^)u + F, (6.67) 

u\0) = S'u"(0) + g. (6.68) 

Since, by assumption, (|6.62p has only the trivial solution for 5ft s > 0, ()6.68p has a unique solution. Inverting the 
Fourier and Laplace transforms gives the solution in physical space. 
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As in the scalar wave case, it is simple to solve (|6.67|) for F = 0. From the inhomogeneous versions of (|6.64p and 

(ma, 

where the CTj are determined by 

C(s, w_)ct = q. 

It is now possible to establish the pseudo-differential version of boundary stability: 

For all cj, s with 3fis > 0, there is a constant K independent of w, s and q such that the solutions of (|6.67p - (|6.68p with 
F = satisfy 

\u{s,{),u)\<K\q{s,u)\. (6.69) 

Boundary stability is equivalent to the requirement that the eigenvalue problem (j6.62p has no eigenvalues for 5ft s > 
or that Det C(aj_, s) 7^ for 5ft s > 0. In particular, it rules out generalized eigenvalues. It is essential to establish the 

Main Theorem: If the half-space problem is boundary stable then it is strongly well-posed in the generalized sense 
of (lOOl) . 

See |46| for the proof, where boundary stability is used to construct a symmetrizer in the Fourier-Laplace representation 
which leads to the estimate 

77||u(s,a;i,cj)f + |u(s,0,l^)P < const (^\\Ff + c\qA . (6.70) 

Inversion of the Fourier-Laplace transform proves the theorem. 

The pseudo-differential theory has far reaching consequences. In particular, the computational rules for pseudo- 
differential operators imply: 

1. The Laplace transform only requires that the estimates hold for 77 > 779 > Oj where 770 is sufficiently large 
to allow for (controlled) exponential growth due to lower order terms. This is essential for extending strong 
well-posedness to systems with variable coefficients. 

2. Boundary stability is also valid if the symbol depends smoothly on (t, x) and is not destroyed by lower order 
terms. Therefore the problem can be localized and well-posedness in general domains can be reduced to the 
study of the Cauchy problem and half-space problems. 

3. The principle of frozen coefficients holds. The properties of the pseudo-differential operators give rise to estimates 
of derivatives in the same way as for standard partial differential equations. Therefore strong well-posedness 
in the generalized sense can be extended to linear problems with variable coefficients and, locally in time, to 
quasilinear problems. 

Since pseudo-differential operators are much more flexible than standard differential operators, they can be applied 
to second order systems as well as first order systems. Consider, for example, the problem (j6.35P " (|6.37p discussed in 
Sec. IVI A 41 After transforming the initial data to zero, the Fourier-Laplace transform becomes 

$,, ^{s'+Lo'')^~F, 
s$ ~a^x - iP^'^ = q- (6-71) 



Introduction of a new variable $a; — -^/jsp + oj^ ^ gives the first order system 



V\ 



1 2, 2 r> ] U — F, U=,-r 

,|2 _J_ 1^2 \ S^ + LJ^ / V * 



* :$~a# /^^ ^ = q, (6.72) 



V|SP+W2 V|s|2+u;2 
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where 

F= , ^ =(l], q= , ^ =q- (6-73) 

The second order problem (|6.35p - (|6.37l) is strongly well-posed in the generalized sense if the corresponding first 
order problem (|6.72p with general data F, q has this property. The second order version of boundary stability ()6.7f|) 
established in Sec. IVI A 4[ rewritten in terms of the first order variables, implies 

ms,0,u:)\ + \i'is,0,u:)\<K\qis,Lo)\, (6.74) 

i.e. the first order order version of boundary stability (|6.69p . Thus the Main Theorem applies and the second order 
problem is strongly well-posed in the generalized sense in the first order version (j6.60p . which is equivalent to the 
second order version (|6.50p . 



2. The energy method for first order symmetric hyperbolic systems 

Energy estimates for first order symmetric hyperbolic systems were first applied to Einstein's equations in harmonic 
coordinates by Fischer and Marsden [53 to give an alternative derivation of the results of Choquet-Bruhat for the 
Cauchy problem. The energy method extends to the quasilinear IBVP with maximally dissipative boundary conditions. 

Again, begin by considering the constant coefficient system ()6.55|) 

m 

ut = V'd.u + F, V'd, ^Ad^.^Y. ^J^-. (6.75) 

on the half-space 

i > 0, xi > 0, — oo < Xj < GO , j ~ 2, . . . ,m, 

with initial data u(0,x) = f{x). As before, A and Bj are symmetric N x N matrices so that the eigenvectors of P* 
form a complete set with real eigenvalues. In matrix notation, there is a symmetric positive definite symmetrizer Hmn 
such that HmpA^ and HmpBj^ are symmetric. Here, the boundary matrix A is allowed to be singular, cf. |49l - l5ll |. 
With an appropriate choice of symmetrizer, it can be put in the form 



-IP 











OQ 











jR 



A^al O^ 0,a>0, u ^ \ u" \ (6.76) 



,u". 



where I^and I^ are identity matrices acting on the P-dimensional subspace u^ and the P-dimensional subspace w^^, 
respectively,and O^ is a zero matrix acting on Q-dimensional kernel u'^ , with N — P + Q + R. No boundary condition 
can be imposed on the components of u^^ , which are the outgoing modes, or the components of u'^ , The components 
of u satisfy PDEs intrinsic to the boundary. They are referred to as zero velocity m,odes since they have no velocity 
relative to the boundary. As in the discussion of the non-characteristic case in Sec. IVIB I[ there are P ingoing modes 
and the IBVP requires P boundary conditions at xi — 0. They are prescribed in the m,aximal form 

u\t,0,x_) = Su^\t,0,x^) + q{t,x^), a;_ = (a;2, . . . ,x„), (6.77) 

where the P x R matrix S satisfies the dissipative condition that, for homogeneous data q — 0, the local energy flux 
out of the boundary be positive, 

T:= {u,Au)>0. (6.78) 

In the simplified form (|6.76p . this leads to the requirement 

- \Su^\t, 0, x-)\^ + \u^^{t, 0, x-)\^ > 0, (6.79) 

where |up = {u, u) in terms of the linear space inner product. 
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The rationale for these maximally dissipative boundary conditions results from the energy estimate for the case of 
homogeneous boundary data q — 0. Beginning with 

m 

dt{u,u) ^2{u,Ad^,u^Y. B,d^^u + F), (6.80) 

integration over the half-space gives 

dtE:^dt\\uf = 2{u,Ad^^u) + 2{u,F) 
= -{u,Au)b + 2{u,F), 
< 2{u,F)<\\uf + \\Ff = E+\\Ff, (6.81) 

where the {u,v) denotes the integral of (m,1'), etc. Thus the maximally dissipative boundary conditions provide the 
required energy estimate. Inhomogeneous boundary data q can be transformed to zero by the change of variable 
M — > u — q to obtain an analogous estimate. 

More generally, if the boundary is uniformly characteristic so that the kernel of the boundary matrix A has constant 
dimension Q, then the quasilinear IBVP problem with maximall y d issi p ative boundary conditions is strongly well- 
posed: a solution exists locally in time which satisfies (I6.60p . See |5ll [SSroq for details. 

The theory can also be recast in the covariant space-time form 

A'^df.u = F (6.82) 

where A'^ — {A*,A^) are symmetric matrices and A* is positive definite. As an illustration of this and the various 
choices of boundary conditions, consider the IBVP for the scalar wave equation 

g'"'d^,d^<i> = 0, x>0, t>0, g""" > 0. (6.83) 

This can be rewritten in the first order symmetric hyperbolic form (|6.82|) for a 5-component field u by introducing 
the auxiliary variables 

M = f Mt 1 = f 5t$ I . (6.84) 

The matrices A'^ are then given by 



A* = ( -.g" I , A' = \ -2.g" -g^' | , (6.85) 

-g'^ 




with 




F^ \ . (6.86) 



In this first order form, the Cauchy data consist of w|t=o = / subject to the constraints 

C^ := u, - a,wo - 0. (6.87) 

The evolution system implies that the constraints satisfy 

dtC, = 0, (6.88) 

so that they propagate up the timelike boundary at x = and present no complication. 

The boundary matrix {—A^) (oriented in the outward normal direction) has a 3-dimensional kernel, whose trans- 
posed basis consists of the zero- velocity modes 

(1,0,0,0,0), (0,0,-g-^g-^0), and (0,0, -5"^0,g--). (6.89) 
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In addition, there is one positive eigenvalue and one negative eigenvalue 

A±-:±A + .g^*, (6.90) 

where 

A = 



^{g-'Y+5,,g-^g-K (6.91) 



Thus precisely one boundary condition is required. 
In terms of the normalized eigenvectors 




e± = -^=== ±\ + g^' . (6.92) 



u — u_|_e+ + u_e_ + uo, where uq lies in the kernel. The boundary condition takes the form u-)_ — Su^ = q, subject 
to the dissipative condition 

J-=-(u,A^u) >0. (6.93) 

For the case of homogeneous data 9 = 0, this requires that 

S^ < -A_ / A+ . (6.94) 



The limiting cases S — ±-y/— A_ / A+ lead to reflecting boundary conditions. In the case of a standard Minkowski 
metric g'^'^ = 77^"^, these correspond to the homogeneous Dirichlet condition and Neumann conditions 

dt<i>\B = 0, d^^B = 0, (6.95) 

and the choice S* = corresponds to the Sommerfeld condition 

(at-c),)$|B = 0. (6.96) 

More generally, the geometric interpretation of these boundary conditions is obscured because the energy E (|6.81|) 
standardly used in the first order theory is constructed with the linear space inner product, as opposed to the 
geometrically defined energy (|6.16p natural to the second order theory. It is possible to reformulate the covariant 
theory of the wave equation in first order symmetric form with boundary conditions based upon the covariant energy 
(|6.16p . But without such a guide to begin with, a first order symmetric hyperbolic formulation can lose touch with 
the underlying geometry. 

This scalar wave example also illustrates how the number of evolution variables and constraints escalate upon 
reduction to first order form. As a result, in the case of Einstein's equations, the advantages of utilizing symmetric 
hyperbolic theory is counterbalanced by the increased algebraic complexity which is further complicated by the wide 
freedom in carry out a first order reduction. See Sec. IVIBI 



C. The characteristic Initial-boundary value problem 

There is another IBVP which gained prominence after Bondi's [3] success in using null hypersurfaces as coordinates 
to describe gravitational waves. In the null-timelike IBVP, data is given on an initial characteristic hypersurface and 
on a timelike worldtube to produce a solution in the exterior of the worldtube. The underlying physical picture is 
that the worldtube data represent the outgoing gravitational radiation emanating from interior matter sources, while 
ingoing radiation incident on the system is represented by the initial null data. 

The characteristic IBVP received little attention before its importance in general relativity was recognized. Ken- 
dall [53] established well-posedness of the double null problem for the quasilinear wave equation, where data is given 
on a pair of intersecting characteristic hypersurfaces. He did not treat the characteristic problem head-on but re- 
duced it to a standard Cauchy problem with data on a spacelike hypersurface passing through the intersection of the 
characteristic hypersurfaces so that well-posedness followed from the result for the Cauchy problem. He extended 
this approach to establish the well-posedness of the double-null formulation of the full Einstein gravitational problem. 
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Kendall's approach cannot be applied to the nuU-timelike problem, even though the double null problem is a limiting 
case. 

The well-posedness of the null-timelike problem for the gravitational case remains an outstanding problem. Only 
recently has it been shown that the quasilinear problem for scalar waves is well-posed [58| . The difficulty unique to 
this problem can be illustrated in terms of the l(spatial)-dimensional wave equation 

{df - dl)<P = 0, (6.97) 

where {t, x) are standard space-time coordinates. The conserved energy 

m -\jd~x {{di<^f + {d,^f^ (6.98) 

leads to the well-posedness of the Cauchy problem. In characteristic coordinates {t — t — x^x — t-\-x), the wave 
equation transforms into 

5t9,$ = 0. (6.99) 

The conserved energy evaluated on the characteristics t = const, 

E[t) = I dx{d^^f, (6.100) 

no longer controls the derivative dt^. 

The usual technique for treating the IBVP is to split the problem into a Cauchy problem and local half-space 
problems and show that these individual problems are well posed. This works for hyperbolic systems based upon 
a spacelike foliation, in which case signals propagate with finite velocity. For (|6.97l) . the solutions to the Cauchy 
problem with compact initial data on i = are square integrable and well-posedness can be established using the L2 
energy norm (|6.98p . 

However, in characteristic coordinates the 1-dimensional wave equation (j6.99p admits signals traveling in the -l-x- 
direction with infinite coordinate velocity. In particular, initial data of compact support $(0, x) = f{x) on the 
characteristic i = admits the solution $ = g{t) + f{x), provided that g{0) — 0. Here g{t) represents the profile of 
a wave which travels from past null infinity (x — > —00) to future null infinity {x -^ -l-oo). Thus, without a boundary 
condition at past null infinity, there is no unique solution and the Cauchy problem is ill posed. Even with the boundary 
condition $(<, —00) = 0, a source of compact support S{t, x) added to (|6.99|) . i.e. 

dtd^^ = 5, (6.101) 

produces waves propagating to a; = +00 so that although the solution is unique it is still not square integrable. 
On the other hand, consider the modified problem obtained by setting $ = e"^^, 

dt{d^ + a)"^ = F , ^{0,x) =€-"'' fix), *(t,-oo)=0, a > 0, (6.102) 

where F = e~°'^S. The solutions to (I6.102p vanish at x = +00 and are square integrable. As a result, the problem 
(I6.102P is well posed with respect to an L2 norm. For the simple case where F = 0, multiplication of (|6.102l) by 
(2a^ -|- dx'^ -\- \dt'^) and integration by parts gives 



^dtjdx({dx^f + 2a'-^A ^^Jdx(2{dt^)dx^-{dt^)A <^Jdx{dx^y 



(6.103) 



The resulting inequality 



for the energy 



dtE < constF (6.104) 



E=^ f dxUdx'i')^ + 2a^^'' 



(6.105) 



provides the estimates for dx^ and ^ which are necessary for well-posedness. Estimates for dt"^, and other higher 
derivatives, follow from applying this approach to the derivatives of (J6.102I) . The approach can be extended to include 
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the source term F and other generic lower differential order terms. This allows well-posedness to be extended to the 
case of variable coefficients and, locally in time, to the quasilinear case. 

Although well-posedness of the problem was established in the modified form (|6.102p , the energy estimates can be 
translated back to the original problem (|6.f0fp . The modification in going from (|6.f0fp to (|6.f02p leads to an effective 
modification of the standard energy for the problem. Rewritten in terms of the original variable $ = e°^^, (|6.105p 
corresponds to the energy 

E=^ f dxe-^""^ ((5.*)' + a'*') • (6.106) 

Thus while the Cauchy problem for (|6.10ip is ill posed with respect to the standard L2 norm it is well posed with 
respect to the exponentially weighted norm (|6.106l) . 

This technique can be applied to a wide range of characteristic problems. In particular, it has been applied to the 
quasilinear wave equation for a scalar field $ in an asymptotically flat curved space background with source S, 

g''^V,V6$ = 5($,5e$,x^), (6.107) 

where the metric g"^ and its associated covariant derivative Va are explicitly prescribed functions of ($,a;"). In 
Bondi-Sachs coordinates [3|, |59| based upon outgoing null hypersurface u — const, the metric has the form 

gf.^dx'^dx'' = -{e^^W - r-^hABW'^W^)du^ ~ 2e^^dudr ~ 2hABW^dudx^ + r^hAsdx^dx^, (6.108) 

where x are angular coordinates, such that {u, x ) — const along the outgoing null rays, and r is an areal radial 
coordinate. Here the metric coefficients {W, /?, W^, Hab) depend smoothly upon ($, u, r, x^) and fall off in the radial 
direction consistent with asymptotic flatness. The null-timelike problem consists of determining $ in the exterior 
region given data on an initial null hypersurface and on an inner timelike worldtube, 

$(0,r,x^) = /(r,a;^), $(u, i?,a;'^) = (j(m,x^), R < r < 00, u>0. (6.109) 

It is shown in |58| : 

The null-timelike IBVP \6.10T^ - ^6.109\] is strongly well-posed subject to a positivity condition that the principal part 
of the wave operator reduces to an elliptic operator in the stationary case. 

The proof is based upon energy estimates obtained in compactifled characteristic coordinates extending to 1+ . 



VII. HISTORICAL DEVELOPMENTS 



Early computational work in general relativity focused on the Cauchy problem. The IBVP only recently received 
serious attention and is still a work in progress. Although there are two formulations where strong well-posedness has 
been established (see See's lIXl and |X|) . several important questions remain. Along the way, there have been partial 
successes based upon ideas of potential value in guiding future work. 



A. The Frittelli-Reula formulation 

The first extensive treatment of the IBVP for Einstein equations was carried out by Stewart [15], motivated at 
the time by the tremendous growth in computing power which made numerical relativity a realistic approach for 
applications to relativistic astrophysics. His primary goal was to investigate how to formulate an IBVP for Einstein's 
equations which would allow unconstrained numerical evolution. Stewart focused upon a formulation of Einstein's 
equations due to Frittelli and Reula [60, ^M j although his approach is sufficiently general to have application to other 
formulations. 

The Frittelli-Reula system was chosen because it is symmetric hyperbolic for certain choices of the free parameters. 
It is based upon the ADM formalism, with metric (j2.13p . in which all second derivatives are eliminated by the 
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introduction of auxiliary variables. There is a 2-parameter freedom in the choice of first order variables. The 
densitized lapse h^'a, where p is an additional adjustable parameter, and the shift are treated as explicitly prescribed 
variables. In addition to the Hamiltonian and momentum constraints ()2.10p . the integrability conditions arising from 
the auxiliary variables lead to 18 additional constraints. Another adjustable parameter controls the freedom of mixing 
the constraints with the evolution system. 

This net result is that the vacuum Einstein equations reduce to a first order evolution system of the form ()6.55p 
consisting of 30 equations governing the metric variables and 22 equations governing the constraints. This is further 
complicated by the lack of geometric or tensorial properties of the evolution variables. Frittelli and Reula analyzed 
the principal part of the system and showed that it is symmetric hyperbolic for certain values of the adjustable 
parameters. 

The well-posedness of the IBVP for such symmetric hyperbolic reductions of Einstein's equations depends upon 
whether proper constraint preserving boundary conditions can be imposed. Stewart analyzed the eigenvalues of the 
boundary matrix for the linearized system using the Fourier-Laplace method described in Sec. lVI A4I For the evolution 
system governing the metric variables, he identified 6 ingoing modes, 6 outgoing modes and 18 zero-velocity modes 
in the kernel which propagate tangential to the boundary. Thus this system requires exactly 6 boundary conditions. 
From the boundary matrix for the system governing constraint evolution, he identified 3 ingoing modes, 3 outgoing 
modes and 16 zero velocity modes, so that 3 boundary conditions are required for constraint enforcement. 

The Fourier-Laplace analysis of the constraint system showed that the determinant condition (|6.66p was satisfied, 
i.e. the homogeneous problem had only the trivial solution, and that the linearized system had a well-posed IBVP. 
The three boundary conditions could be satisfied by requiring that the Cauchy momentum constraints (j2.17p vanish 
on the boundary. 

The analysis of the evolution system governing the metric showed that the constraints could be enforced by a 
particular choice of boundary data for the metric variables. In this way, the evolution could be freed from the 
constraint system. The details are hidden in the symbolic algebra scripts necessary to analyze the complexity of 
the Fourier-Laplace modes. No discussion was given of the estimates for the derivatives which would be necessary 
for the well-posedness of the nonlinear problem. There has apparently been no further results for the Frittelli-Reula 
formulation and no attempt at a numerical evolution code. 

Although the results were inconclusive for the nonlinear case, Stewart's treatment provided the first example of 
how to apply pseudo-differential techniques to the IBVP for Einstein's equations and served as the basis for much 
of the following work. The recognition that constraint preservation for this system could be achieved by enforcing 
the Cauchy momentum constraints on the boundary suggests a possible wider application but whether there is any 
universal procedure for enforcing the constraints in the 3 + 1 formulation remains an open issue. See Sec. IXII for a 
discussion. 



B. The BSSN and NOR formulations 



Codes based upon the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation [g^, |6j| have been used by the 
majority of groups [6J-|68| carrying out simulations of binary black hole and neutron star systems. The successful 
development of the BSSN formulation proceeded through an interplay between educated guesses and feedback from 
code performance. Only in hindsight has its success spurred mathematical analysis, which showed that certain 
versions were strongly hyperbolic and thus had a well-posed Cauchy problem [69-71J. Although significant progress 
has been made in establishing some of the necessary conditions for well-posedness and constraint preservation of the 
IBVP [72, U3, UM 1 there is still no satisfactory mathematical theory on which to base a numerical treatment of the 
boundary. 

Similar to the Frittelli-Reula system, the BBSN system reduces the Einstein equations to first order form by the 
introduction of auxiliary variables. In addition, there is a conformal-traceless decomposition of the 3-geometry. The 
Nagy-Ortiz-Reula (NOR) system |70| is similar but without the conformal decomposition. Again, a number of free 
parameters enter into the first order reduction and into the way that the constraints are mixed with the evolution 
system. For a particular choice of parameters the linearization off Minkowski space is symmetric hyperbolic [75| and 
leads to a well-posed IBVP for the linearized problem. However, the corresponding nonlinear system is no longer 
symmetric hyperbolic and there is no well-posed boundary treatment. 

Gundlach and Martin-Garcia |72l] studied simplified second order versions of of the BSSN and NOR systems which 
were symmetric hyp erbolic in a sense they defined in [73|. They were able to confirm and generalize many of the 
results found in [63| for the first order reduction. Of particular interest, they found that all the characteristic modes 
propagated causally, in contrast to the superluminal modes present in the first order system. The chief shortcoming 
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of their treatment is the incompatibihty of the constraints with the dissipative boundary conditions necessary for 
well-posedness. 

The strong well-posedness of the IBVP for 3 + 1 formulations remains an outstanding problem. The strategy 
in current numerical practice for BSSN evolution systems is to apply naive, homogeneous Sommerfeld boundary 
conditions, where needed, to each evolution variable (cf. [63]) and place the boundary out far enough so that its 
harmful effects are limited. 



Other 3 + 1 studies 



Many of the other early investigations on the well-posedness of the IBVP for 3 + 1 formulations centered about 
linearized systems |77H79l l83|. spherically symmetric and ID spacetimes [84l - [87| or other model problems |73l. ISSl - IOOj 
which simplified the treatment. In particular, well-posedness of the linearized problem is a necessary condition for 
extension to the nonlinear case. 

One promising approach was based upon a generalization of the Frittclli-Reula and Einstein- Christoffel (EC) [91[ 
systems, which Kidder, Scheel and Teukolsky (KST) showed was symmetric hyperbolic for certain values of the free 
parameters |92| . In [78], energy estimates for maximally dissipative boundary conditions were used to formulate a 
well posed IBVP for the linearization of this system off Minkowski space. However, constraint preservation limited 
the allowed boundary conditions to the reflecting Dirichlet or Neumann type. 

The geometric analogy between the Hamiltonian and momentum constraints (|2.10p of the Cauchy problem and the 
boundary constraints 

G°''iVb = 0, (7.1) 

where N}, is the normal to the boundary, led Frittelli and Gomez to propose that (|7.ip be enforced as boundary 
conditions 804821 • They showed for the EC system, with vanishing shift and certain choices of the free parameters, 
that enforcing three linearly dependent combinations of the boundary constraints would lead to preservation of the 
Hamiltonian and momentum constraints. Furthermore, they showed that these linear combinations could be used to 
formulate boundary conditions for three of the ingoing metric variables. They did not study the boundary stability 
of the resulting IBVP. 

In [831 , the determinant condition (I6.66[) of the Fourier-Laplace method was applied to the linearized (EC) system to 
identify ill-posed modes arising from various choices of boundary conditions and free parameters. Several noteworthy 
results were found. For parameter choices in which the EC system was strongly but not symmetric hyperbolic, they 
found that maximally dissipative boundary conditions gave rise to ill posed modes. This is in accord with the general 
theory which requires both maximally dissipative boundary conditions and a symmetric system to guarantee a well 
posed IBVP. In addition, for a range of parameters giving rise to a symmetric system, ill-posed modes were found for 
boundary conditions based upon the boundary constraints (I7.ip . As in the case of the Cauchy momentum constraints 
proposed by Stewart, this casts doubt on whether such boundary constraints are universally applicable. The approach 
in [83] was effective for ferreting out what doesn't work but did not go beyond the results of [73] for establishing a 
well-posed IBVP. 

In a later study [93[, the EC system was further generalized to include a dynamical lapse of the Bona-Masso 
type [9J] and fixed shift. The IBVP was analyzed in the high frequency limit, again using the determinant condi- 
tion of the Fourier-Laplace method to determine ill-posed modes. It was found that constraint preserving boundary 
conditions that were based upon the Newman-Penrose [9a] Weyl curvature component ^o satisfied the determinant 
condition provided the evolution system was strongly hyperbolic and the constraint propagation system was sym- 
metric hyperbolic. Boundary conditions based upon the Newman-Penrose ^o component were first introduced in the 
Friedrich-Nagy system [l^. Other ^q boundary conditions for the EC system were tested in |79l. l89l. 1961- 100. 102]. 
They have been improved to be highly effective absorbing boundary conditions for gravitational waves (see Sec. IVIIip 
and have performed well in numerical tests. See Sec. IXIIIl 



D. The harmonic and Z4 formulations 



The IBVP for general relativity takes on one of its simplest forms in the harmonic formulation, in which the Einstein 
equations reduce to 10 quasilinear wave equations. Nevertheless, progress on this problem was not straightforward. 
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Difficulties arose in liandling tlie fiarmonic constraints (I2.2p . in wliicfi derivatives of tfie metric tangential to the 
boundary prohibited use of standard dissipative boundary conditions for the wave equation. An e arly well-posed 
treatment was based upon the observation that the harmonic Cauchy problem is well-posed [ij, Il03j | . Consequently, 
if locally smooth reflection symmetry were imposed across the boundary then the well-posedness of the Cauchy 
problem would imply well-posedness of resulting IBVP on either side of the boundary. The reflection symmetry forces 
the troublesome tangential derivatives to vanish but it also forces homogeneous boundary conditions of the Dirichlet 
or Neumann type. Although these boundary conditions satisfy the dissipative criterion for a well-posed IBVP, they 
were too restrictive for use in practical numerical applications and did not allow large boundary data. It took a 
different approach to formulate a strongly well-posed harmonic IBVP with Sommerfeld type boundary conditions. 
See Sec. H 

The ZA, formalism [104| aims at a covariant version of hyperbolic reduction by expressing the vacuum Einstein 
equations in the form 

QI.V _ y(M^^) + i^M^v^^P = 0. (7.2) 

When the vector field Z^ is identified with the (generalized) harmonic conditions C^, this reduces to the harmonic 
formulation. However, the freedom is retained to introduce other gauge co nditions which force Z^ — 0. When only 
6 components of (|7.2|) are required to vanish, it has been shown [271 . llOSj that the ^4 formalism encompasses the 
standard 3 -f 1 formulations, including the ADM, NOR, BSSN and KST systems. 

It is possible that the close analogue between the Z4 and harmonic formulations might be used to shed li ght on th e 
IBVP for 3-1-1 systems. Constraint preserving boundary conditions for such Z4 systems have been proposed 106l4l08l |. 



However, as for other 3 + 1 formulations, the boundary stability necessary for a strongly well-posed IBVP has not 
been established. Nevertheless, the results of numerical tests are promising. See Sec. IXIIII 



VIII. NON-REFLECTING OUTER BOUNDARY CONDITIONS 

The correct physical description of an isolated system involves asymptotic conditions at infinity which ensure that 
the radiation fields have the proper 1/r falloff and that the total energy and radiative energy loss are finite. This can 
be achieved by locating the outer boundary at T^ . Instead, current simulations of binary black holes are carried out 
with an outer boundary at a large but finite distance in the wave zone, i.e. many wavelengths from the source. This 
is in accord with the standard practice in computational physics to impose an artificial boundary condition (ABC) 
which attempts to approximate the proper behavior of the exterior region. 

At the analytic level, many ABCs are possible, even Dirichlet or Neumann conditions, provided the proper boundary 
data is known to allow outgoing radiation to pass through. However, the determination of the correct boundary data 
is a global problem, which requires extending the solution to X+ either by matching to an exterior (linearized or 



nonlinear) solution obtained by some other means. The matching approach has been reviewed elsewhere IJ]. As 
shown by Gustafsson and Kreiss, the construction of a non-reflecting b ound ary condition for an isolated system in 
general requires knowledge of the solution in a neighborhood of infinity |109| . 

Even if the outgoing radiation data for the analytic problem were known, at the numerical level a Dirichlet or 
Neumann condition would reflect waves generated by the numerical error and trap them in the grid. The alternative 
approach is an ABC which is non-reflecting for homogeneous data. Artificial boundary conditions for an isolated 
radiating system for which hom ogeneous data is approximately valid are commonl y ca, l led a bsorbing boundary con- 
ditions (see e.g. IllOL Ill2l4ll6ll '). or non-refiecting boundary conditions (see e.g. |ll7l - ll2Clj |) or radiation boundary 



conditions (see e.g. (l24 Il25| ). Such boundary conditions are advantageous for computational use. However, local 
ABCs are in general not p erfect. Typically they cause some partial reflection of an outgoing wave back into the 
system [112| . Ill6l Il26l Il27j . It is only required that there be no spurious reflection in the limit that the boundary 
approaches an infinite sphere. 

A traditional ABC for the wave equation is the Sommerfeld condition. For a scalar field $ satisfying the Minkowski 
space wave equation (|6.ip with compact source, the exterior retarded field has the form 

9 — \ ! , (S.ij 

where /, g and h are smooth bounded functions. The simplest case is the monopole radiation field 

$ = l^tlA (8.2) 



31 

which satisfies {dt + dr){r^) — 0. This motivates the use of the Sommerfeld condition 

-{dt+dr)ir'^)\R^q{t,R,0,(b) (8.3) 

r 

on a finite boundary r = R. However a homogeneous Sommerfeld condition, i.e. g = 0, is exact only for the spherically 
symmetric monopole field. The Sommerfeld boundary data q for the general case (|8.1|) falls off as 1/R'^, so that a 
homogeneous Sommerfeld condition introduces an error which is small only for large R. As an example, 



/(i-i?)cos6' 
<i = 

for the dipole solution 



»= " a3 (8") 



*„„.,. 4Z<i^.-f£<i^ + /<i^')c„.<,, (8^5) 



A homogeneous Sommerfeld condition at r = i? leads to a solution ^Dipoie containing a reflected ingoing wave. For 
large R, it is given by 

;;. ^ F{t + r ~2R)cose .^ . 

^Dipolc - *Dipolc + K^ , (8.6) 

where dtf{t) = F{t) and the reflection coefficient has asymptotic behavior k = 0{1/R^). More precisely, the Fourier 
mode 



$Dipoic(w) = aj + K^ j (8.7) 

satisfies the homogeneous boundary condition [dt + dr){'''^ Dipoie{^))\R — with reflection coefficient 

11 



k(cj) = 



2a;2i?2 + 2iujR - 1 2lj^R'^ ' ' ' ' 

Note, from dHU) and dH]), 

K ~ qR. (8.9) 

Use of this relationship simplifies the determination of the asymptotic falloff of the reflection coefficient by avoiding an 
explicit calculation of the reflected wave. Also note, that if (|8.3p is replaced by the second order Sommerfeld condition 

\{dt + dry{dt + dr){r<P)\R = 92 (8.10) 

then dipole radiation has homogeneous data q2 — 0. In this way, the falloff rate of the reflection coefficient is reduced 
from dSS]) to K2 ~ l/R^- 

The exponent n of the 0(l/i?") falloff of the refiection coefficient is an important measure of the accuracy of an 
ABC. Such reflection coefficients can be calculated for linearized gravitational wav es on a M inkowski background, 
analogous to the above scalar wave example, using either the R egge-Wheeler-Zerilli |128l4l30| perturbative method, 
as carried out in [99, 100 , 102] , or by the Bergmann- Sachs |l3l| gravitational Hertz potential method, as carried out 
in [i3|. See Sec. IX Al The main difference in the gravitational case arises from the gauge modes, which exist along 
with the radiative degrees of freedom. In first order formulations, this is further complicated by the modes introduced 
by the auxiliary variables. The second order harmonic formulation, in which all modes propagate on the light cone, 
is simplest to analyze. See Sec. |X] for a discussion of reflection coefficients in the h armonic case. 

Local ABCs have been extensively applied to linear problems with varying success J110l . lll2| - |ll5lll24ll26j |. Son ie ar e 
local approximations to exact integral representations of the solution in the exterior of the computational do main 
while others are based on approximating the dispersion relation of the so-called one-way wave equations 113l 



Higdon [112l | showed that this last approach is essentially equivalent to specifying a finite number of angles of incidence 
for which the ABCs yield perfect transmission. Local ABCs hav e also been derived for the linear wave equation by 
considering the asymptotic behavior of outgoing wave solutions 'l24!|, thus generalizing the Sommerfeld condition. 
Although this type of ABC is relatively simple to implement and has a low computational cost, the final accuracy can 
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be limited if the assumptions about the behavior of the waves are oversimphfied. See [ll6l IllSl . Il32l Il33j | for general 
discussions. 

The disadvantages of local A BCs have led t o co nsideration of nonlocal versions based on integral representations 
of the infinite domain problem 118 yi20l Il32l Il34| . Even when the Green function is known, such approaches were 
initially dismissed as impractical iHOJ; however, the rapid development of computer power and numerical techniques 
has made it possible to implement exact nonlocal ABCs for the linear wave equation and Maxwell's equations in 
3D [l35|, |l36|. If properly implemented, this method can yield numerical solutions to a linear problem which converge 
to the exact infinite domain problem in the continuum limit, while keeping the artificial boundary at a fixed distance. 
However, due to nonlocality, the computational cost per time step usually grows at a hi gher power with grid size, 
0{N'^) per time step in three spatial dimensions, than for a local approach J118l . ll32ill35j |. 

The extension of ABCs to nonlinear problems is much more difficult. The problem is normall y treated by linearizing 
the region between the outer boundary and infinity, using either local or nonlocal linear ABCs (l32lll33| . The neglect 
of the nonlinear terms in this region introduces an unavoidable error at the analytic level. However, even larger errors 
are typically introduced in prescribing the outer boundary data. The correct boundary data must correspond to the 
continuity of fields and their normal derivatives when extended across the boundary into the linearized exterior. This 
is a subtle global requirement for any consistent boundary algorithm, since discontinuities in the field or its derivatives 
would otherwise act as a spurious sheet source on the boundary that would contaminate both the interior and the 
exterior solutions. However, the fields and their normal derivatives constitute an over determined set of data for the 
boundary problem. So it is necessary to solve a global linearized problem, not just an exterior one, in orde r to find 
the proper data. An expedient numerical method to eliminate back reflection is the use of sponge layers, cf. |12l| , in 
which damping terms are introduced into the evolution equations near the outer boundary. 

The designation "exact ABC" is given to an ABC for a nonlinear system whose only error is due to linearization of 
the exterior. An exact ABC requi res the us e of global techniques, such as the boundary potential method, to eliminate 
back reflection at the boundary |120l Il32l |. Furthermore, nonlinear waves intrinsically backscatter, which makes it 
incorrect to try to e ntirely elim inate incoming radiation from the outer region. For the nonlinear wave equation, test 
results presented in |l22l Il23j | showed that Cauchy-characteristic matching outperformed all ABC's in the existent 
literature. 

It is an extra challenge to apply ABCs to strong ly n onlinear hydrodynamic problems [118| . Thompson |l37 



generalized a previous nonlinear ABC of Hedstrom [117{ to treat ID and 2D problems in gas dynamics. These 



boundary conditions performed po orly in so me situations because of d ifficu lties in adequately modeling the field 
outside the computational domain [llSl Il37| . Hagstrom and Hariharan [l38J | have overcome these difficulties in ID 
gas dynamics by their use of Riemann invariants. They proposed, at the heuristic level, a generalization of their local 
ABC to 3D. 

In order to reduce the analytic error, an artificial boundary for a nonlinear problem must be placed sufficiently far 
from the strong-field region. This can increase the computational cost in multi-dimensional simulations [110]. There is 
no ABC which converges (as the discretization is refined) to the infinite domain exact solution of a strongly nonlinear 
wave problem in multi-dimensions, while keeping the artificial boundary fixed. When the system is nonlinear and not 
amenable to an exact solution, a finite outer boundary condition must necessarily introduce spurious effects. Attempts 
to use compactified Cauchy hypersurfaces which extend the dom ain t o spatial infinity have failed because the phase 
of short wavelength radiation varies rapidly in spatial directions [127]. In fact, in his pioneering simulation of binary 
black holes, Pretorius [28, 29] used this effect as a numerical expedience by applying artificial dissipation to diminish 
short wavelength error arising from the use of a compactifi ed o uter boundary at spatial infinity. For a recent review 
of ABCs in the computational mathematics literature, see (l39J | . 

The situation for the gravitational IBVP is not as severe as for hydrodynamics, especially for formulations in which 
the gauge modes and radiation modes propagate with the same speed. However, due to nonlinearities, there is always 
some error of an analytic nature introduced by a finite boundary which is independent of discretization. In general, a 
systematic reduction of this error can only be achieved by moving the computational boundary to larger and larger 
radii. There has been recent progress in designing absorbing boundary conditions for the gravitational field. Buchman 
and Sarbach l99| have developed higher order local boundary conditions based upon derivatives of $o which are non- 
reflecting up to any given multipole mode for linearized gravitational waves, analogous to the switch from a first order 
Sommerfcld condition (|8.3p to the second order condition (|8.10p . These boundary conditions freeze the value of ^o 
to its initial value, i.e. 

5**0 = 0, (8.11) 

in order to avoid a Oth order violation of the compatibility condition between the initial and boundary data. They 
have extended this approach to include quadrupole gravit ation al waves on a Schwarzschild backg round and also to 
account for 0{M/R) back scatter using a nonlocal version [lOOl ]. For a review of this approach see [lOll ]. It has been 
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applied to the harmonic formulation in (l02J | and to the Z4 formulation in [l08J | . These are possibly the best possible 
local b oundary c onditions for the gravitational radiation modes. 

Lau 14Cll - ll42J has formulated a n ex act ABC for linearized gravitational waves on a Schwarzschild background. 



Based upon the flat space work of |12Cll | , he reduces the calculation of the Green function incorporating the boundary 
condition for the perturbed metric to the integration of a radial ODE by using a combined spherical harmonic and 
Laplace transform of the Regge- Wheeler- Zerilli equations. He discusses the trade-off in computational cost for this 
nonlocal ABC versus the larger computational domain required by a local condition. This is similar to the trade-off 
between characteristic matching and the application of local boundary conditions. 



IX. THE FRIEDRICH-NAGY SYSTEM 



Friedrich and Nagy [lal have presented a theorem establishing the first strongly well-posed IBVP for Einstein's 
equations with the generality to handle an outgoing radiation boundary condition. The approach uses the energy 
method for first order symmetric hyperbolic systems described in Sec. IVIB21 Their formulation is based upon the 
Einstein-Bianchi system of equations with evolution variables consisting of an orthonormal tetrad ea, a = (0, 1, 2, 3), 
the associated connection coefficients F? and Weyl curvature tetrad components C-ggj. Although it differs from 
the metric based formulations used in numerical relativity, the success of their treatment suggests that many of the 
underlying ideas should be universally applicable. In their words, "There are certainly many possibilities to discuss 
the initial boundary value problem and there will be as many ways of stating boundary conditions. However, all of 
these should be just modifications of the boundary conditions given in our theorem" . Perhaps this is overstated since 
their formulation is 3rd differential order in the metric, as opposed to the 2nd order 3 -I- 1 or harmonic formulations. 
Yet, all successful formulations must have the common property of prescribing data which produces a unique solution 
to Einstein's equations. 

The tetrad vector ep is chosen to be timelike and tangent to the boundary. It is used to construct adapted 
coordinates x'^ — {t,x^) = {t,x,y,z) satisfying 

C^,t = 1 , C^,x' = 0, (9.1) 

so that eg plays the role of the evolution vector field t^ introduced in Sec. IIVI However, the evolution field now has 
the metric property of being a timelike unit vector so, as a reminder, I denote it by T'^ = Cg. In accord with the 
notation in Sec. IIV[ let the initial hypersurface So be given by i = and the boundary T be given by a; = 0, with 
adapted coordinates {t,x^). The tetrad vectors are adapted to the geometry as follows. Let N^^ — ei^ oc — V^o; be 
the unit outer normal to T. 

Extend N'^ throughout the spacetime manifold Ai by requiring that it be the unit normal to the hypersurfaces Tc 
given hy X = c = const > 0. On Sq, the remaining tetrad vectors e^, A = (2,3), are chosen to be an orthonormal 
dyad for the {t — 0,x = c) subspaces. They are then propagated throughout A4 by Fermi- Walker transport along 
the integral curves of T'^, which lie in Tc- The connection components intrinsic to Tc are considered to be freely 
specifiable gauge source functions. In addition, the mean extrinsic curvature K oiTc is also considered to be a gauge 
source function. Moreover, it is shown that the equation K — q{x^) can be cast as a quasilinear wave equation which 
determines Tc given initial data corresponding to a; = c and dtx = at i = 0. 

By design, this choice of adapted coordinates and tetrad gauge greatly simplify the IBVP. The evolution system 
consists of the gauge conditions governing the tetrad, the equations relating the connection to the tetrad, the vacuum 
equations relating the Weyl curvature to the connection (which imply the vanishing of the Einstein tensor) and the 
vacuum Bianchi identities, i.e. the tetrad version of the equations 

V^C^.p, = 0. (9.2) 

The system is over determined and subject to constraints arising from integrability conditions. Remarkably, it can 
be reduced to a system with the properties: 

• The evolution system is symmetric hyperbolic. 

• The boundary matrix (j6.76p admits two ingoing variables corresponding to combinations of the vI/q and ^4 
Newman-Penrose components of the Weyl tensor. 
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• Maximally dissipative boundary conditions take the form 

*o + a*4 + /3*4 = g, (9.3) 

where q is the (complex) boundary data and a and /3 are coefficients subject to a dissipative condition. (Con- 
ventions here are chosen to be consistent with Newman-Penrose conventions and differ from 16].) 

• The subsidiary system governing the constraints is symmetric hyperbolic and intrinsic to the Tc hypersurfaces, 
i.e. all derivatives are tangential to Tc- This gives rise to constraint preservation without requiring any further 
restrictions on the boundary conditions. 

The resulting IBVP is strongly well-posed. Given initial data on S, including the hyperbolic angle Q in (|3.4[) at 
the edge Bq, the boundary data q, a choice of gauge source functions and a dissipative choice of boundary condition, 
there exists a unique solution locally in time. Furthermore, the solution depends continuously on the data. 

Several important points should be noted: 

• The specification of the mean extrinsic curvature K of the boundary geometrically determines the location of 
the boundary. 

• The choice of unit timelike vector T" tangent to the boundary represents gauge freedom in the construction of 
the solution. It induces a corresponding foliation of the spacetime and the boundary according to ^C^t — 1. 

• The geodesic curvature of the integral curves of T"* constitute gauge source functions required on the boundary. 
This gauge freedom feeds into the adapted coordinates (i,x^) of the boundary. As a result, the functional 
specification of K{t,x^) becomes gauge dependent. This complication could be avoided by choosing T°- to 
be geodesic but at the expense of possible coordinate focusing singularities which would affect the long term 
existence of the solution. 

• The outgoing null vector K"^ and ingoing null vector L" used in defining ^I^o and ^1^4 are determined by the choice 
of T° (gauge) and the boundary normal TV" by 

K'^^T'^ + N'', L'' = T''-N''. (9.4) 

Friedrich and Nagy are careful to point out that gauge freedom prevents any meaningful interpretation of ^I'o 
and ^i as either purely outgoing or ingoing radiation. 



X. THE HARMONIC IBVP 



The first successful treatment of the harmonic IBVP for Einstein's equations was carried out using the pseudo- 
differential theory described in Sec. I VI A4l which established strong well-posedness in a generalized sense [l7|. The 
theory was developed for the second order formulation of the generalized harmonic formulation (|2.2I) - (|2.5I) . The 
boundary conditions were given in Sommcrfeld form in terms of the outgoing null vector K"" — T" -I- N"" normal to 
the foliation of the boundary. Here, as in the Friedrich-Nagy approach (|9.4I) . T" is a future directed timelike unit 
vector tangent to the boundary but now it is also chosen to be normal to its foliation Bt- Recall that n°, the normal 
to the Cauchy foliation iSt, is not necessarily tangent to the boundary. The motion of the boundary, characterized by 
the hyperbolic angle p.4p . distinguishes T° from n°. 

In the adapted coordinates x^ — {t>Q^x> Q^x^) described in (14.31) . six Sommerfeld boundary conditions for the 
densitized metric 7^"^ = ^/—gg^'^ were given by 

i^^a^7^^ = q^''{t,x'') (10.1) 

i^'^5^(7*^ + 7"^) - q^{t.x^) (10.2) 

if'^947" + 27*"+7"") = q{t.x^), (10.3) 
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where the q's are freely prescribed Sommerfeld data. The harmonic constraints (|2.2I) were used to supply four 
additional boundary conditions, which could be expressed in the Sommerfeld form 

+ ^(9t + 9.)(7*^+7"'^) + aB7^''-\/^r^-0 (10.4) 

V^{C'+n = 1(9, ^ 9,) (f*-^-) 

+ ^{dt + 9x)(7" + 27"* + in + Ml'^ + 7"^) - A/^(f * + f-) = (10.5) 

V^C = i(a,-9,)(7"-7*") 

+ \{dt + 9x)) Ul" + 27*" + in + (7" - 7"")) + dAl*^ - y^f * = 0. (10.6) 

Constraint preservation then follows from the homogeneous wave equation p.6p . 

The key feature of (I10.ip - (|10.6p is that they form a sequential hierarchy of Sommerfeld boundary conditions for the 
metric variables such that the source terms are given in terms of derivatives of previous variables in the sequence. For 
instance, the terms on the second line of (|10.4p are derivatives of 7"^^ and (7*"^ + 1^^), whose boundary conditions 
are prescribed previously in (jlO.ip and ()10.2p . This pattern persists for the remaining boundary conditions in the 
sequence, i.e. (|10.5I) and (|10.6I) . This structure gives rise to a corresponding sequence of estimates for the variables in 
the hierarchy, which is the key for establishing boundary stability and the strong well-posedness of the IBVP. There 
is considerable freedom in the boundary conditions provided this hierarchical structure is preserved. 

The well-posedness of the harmonic IBVP was subsequently also established using estimates for the non-standard 
energy (I6.16P associated with a timelike vector pointing outward from the boundary !18|. The hierarchical structure 
of the boundary conditions corresponds to the upper triangular property (|VI A 31) , which is sufficient condition for a 
well-posed IBVP for a coupled system of quasilinear wave equations. _ 

A more general and geometrical version of these results in terms of a background metric g^^^ was presented in [19| . 
The connection Vo and curvature tensor R'^cab associated with the background metric g^^^, have the same tensorial 
properties as the corresponding quantities Vq and R'^cab associated with gab- In particular, the difference Vq — Vi 
defines a tensor field Cf^ according to 

(Va - Ky = c^y (10.7) 

for any vector field v''. In terms of the (nonlinear) perturbation 

fab = 9ab - 9ab (10-8) 

of the metric from the background. 



c 



ab 



^.9^' (4/6c + Va/bc - yjab) ■ (10.9) 

Since (?^{, is explicitly known, a solution for fab is equivalent to a solution for gab- Einstein's equations are given by 

^ab ._ Qab _ ^(aQb) ^ }_ ^ab^ ^Qd ^ q (10.10) 

subject to the harmonic constraints 

C ■■= g'^^'Ct, = 0. (10.11) 
In the adapted coordinates, the harmonic constraints take the form 

C'':=.g^'^(r^.-r^.) = o, (10.12) 

so that the background Christoffel symbols T^ appear as harmonic gauge source functions. When the harmonic 
constraints are satisfied, the reduced Einstein equations form the desired quasilinear wave system for f^y, 

9'"'yp^af^.. = 2gxrg'"'C\pC\, + W,^^gy)xC\r9''^ - 2g''"i?\„(^5,)A. (10.13) 
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The analogue of the Sommerfeld conditions (|10.ip - (|10.6p are prescribed in terms of the boundary decomposition of 
the metric 

gab = NaNb + Hab , Hab = ^TaTb + Qab, (10.14) 

This leads to an orthonormal tetrad (T", A^°, Q", Q") on T, where Q" is a complex null vector tangent to Bt with 
normalization 

Qab^QiaQb), Q"Qa = 2, Q''Qa=0. (10.15) 

(The tetrad is unique up to the spin freedom Q" — > e*^Q° which does not enter in any essential way.) In terms of the 
outgoing and ingoing null vector fields K"^ — T'^ + N"' and L"^ — T"" — N"-, respectively, which are normal to Bt, the 
metric has the null tetrad decomposition 

gab = -Kf^aLb) + Q(aQb)- (10.16) 

Six Sommerfeld boundary conditions which determine the components of the outgoing null derivatives K'^Vafbc are 
then given by 

Ix'-K-K'^Vahc = q'^Ka, (10.17) 

{Q^'K'K^ - ^K''K^Q^)\/afbc - g"Qa, (10.18) 

{L'K^'K'' - ^K'K'L-)Vafbc = q'^La, (10.19) 

{\Q''Q'K'^~Q''K'Q'')Vafbc = 2(7, (10.20) 

in terms of boundary data q°' and a. The harmonic constraints provide four additional boundary conditions which, 
in terms of the null tetrad, can be expressed in the Sommerfeld form 

- 2C°'Ka = {Q^Q'^K'' + K^K'^L'' - K^Q'^Q" - K^'Q^Q") Vahc = 0, (10.21) 

-2C''Qa = {L^'Q^K'^ + K^'Q^L'^ ~ K^L'^Q'' + Q'-QV") ^afbc = 0, (10.22) 

-2CU, = {L^'L^K" + Q^'Q^L'^ - Q^'L'Q^ - Q'^L^Q'') Vafbc - 0. (10.23) 



As before, constraint preservation follows from p.6p . 

Together, (|10.17p - (|10.23p provide Sommerfeld boundary conditions for the components of K'^Vafbc in the sequential 
order (KK), (QK), (LK), (QQ), (QQ), (LQ), (LL) in terms of the boundary data and the derivatives of preceding 
components in the sequence. This hierarchy of Sommerfeld boundary conditions satisfies the requirements given in 
Sec. IVI A3I (Theorem 1 of [1^) for a strongly well-posed IBVP for the quasilinear hyperbolic system (jlO.131) . See 
Sec. lXIII for a geometrical interpretation of the boundary data and Sec. lXIlIl for numerical tests. 



A. Application to an isolated system 



The main application of the gravitational IBVP is to the spherical outer boundary used in the simulation of an 
isolated system emitting radiation. As discussed in Sec. IVIIH in the absence of an exterior solution, the simplest 
approach is the use of a Sommerfeld boundary conditions with homogeneous data. In doing so it is important to take 
advantage of the freedom in the form of the boundary conditions in order to reduce back reflection. 

Sommerfeld boundary conditions consistent with a well posed harmonic IBVP have wide freedom regarding the 
addition of (i) partial derivative terms consistent with the hierarchical structure and (ii) lower order algebraic terms. 
Various choices were considered in [13| . They were tested by computing the resulting reflection coefficients for spherical 
waves in a Minkowski space background. For this purpose the densitized metric is approximated to linearized accuracy 

by 



-ggf'" ^rf +^t'''^ (10.24) 
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where rj'^'^ is the Minkowski metric. The calculation of the reflection coefficients proceeds as for the scalar wave 
example in Sec. IVIIIi as modified to deal with the gauge modes. 

Li near ized waves in the harmonic gauge can be constructed from the gravitational analogue of the Hertz poten- 
tial [l3l| . which has the symmetries 



and satisfies the flat space wave equation d'^daH'^""^ = 0. Then the perturbation 

satisfies the linearized Einstein equations d^dcr^^" = in the harmonic gauge d^^^" — 0. Outgoing waves can be 
generated from the potential 

^M"^/3 ^ j^p-cufs fii-r) ^ ^^^ ^ i^A'"^/35^5 /(i:ili^ (10.25) 

r r 

where K^^""^^ is a constant tensor. (All higher multipoles can be constructed by taking spatial derivatives.) K^°"^^ has 
2f independent components but the choice K^^""^^ — ^t^°"^P leads to 7^"^ = so there are only 20 independent waves. 
These can reduced to f pure gauge waves for which the linearized Riemann tensor vanishes, which correspond to the 
trace terms in K^°"^l^; e.g. K'^"'^^ = rf'^r]^'^ — rj'^'^rj"'^ leads to a monopole gauge wave. The trace-free part gives rise 
to 10 independent quadrupole gravitational waves, corresponding to spherical harmonics with {(. = 2, —2 < ttt, < 2) in 
the two independent polarization states. 

For a boundary at r = i?, the Sommerfeld derivative in the outgoing null direction is 

K^'d^.^dt + dr. (10.26) 

In formulating boundary conditions which minimize back reflection, the property K^d^j,f{t—r) = is used to introduce 
the appropriate powers of r, analogous to the scalar example (J8.3I) . In [19[, the optimal choice was found to be the 
Sommerfeld hierarchy 






■^K^KpK'^d^ir^j''^) = QKK , (10.27) 



^K^QfjK^d^ir'j'^^) = QKQ , (10.28) 



V 



-Q^Q^.K'^d^ir'r^)-! = Qqq, (10.29) 



r 



QaQpK'^df.j''^ - QaKpQf'df.j''^ = qQQ . (10.30) 

It was found that the data g.. = 0{1/R'*') for the outgoing gravitational quadrupole waves and q,, = 0{\/B?) for the 
outgoing gauge waves. This implies, in accord with (|8.9p . that homogeneous Sommerfeld data gives rise to reflection 
coefficients k = 0{1/R^) for the gravitational waves and k — 0{1/B?) for the gauge waves. These results were 
confirmed us ing the Regge- Wheeler- Zerilli perturbative formulation along with the metric reconstruction method 
described in jl30| |. 

The analysis of linearized waves shows that qqq controls the amplitude of the gravitational radiation passing through 
the boundary. Higher order boundary conditions can be based upon replacing (|10.30p by a condition on ^I^o, which 
in the linearized theory controls the radiation in a gauge independent manner. In this way, the V&o based boundary 
conditions discussed in Sec. lVIlil can be used to further increase the 1/i?" falloff rate of the reflection coefficients for 
gravitational waves. 



XI. CONSTRAINT PRESERVATION 



The IBVP for Einstein's equations is still not well understood due to a great extent from complications arising from 
the constraints. The Hamiltonian and momentum constraints on the Cauchy data take the universal form (j2.10p in 
terms of the components of the Einstein tensor normal to the initial hypersurface. However, there is no common way 
to ensure constraint preservation for the various formulations of Einstein's equations. Even the constraints themselves 
take on different forms. 
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The Friedrich-Nagy system (see Sec. IIX[) is based upon the Einstein-Bianchi equations which is third differential 
order in terms of the metric. In that case, by a cleverly designed choice of adapted coordinates and gauge, constraint 
propagation is governed by a system tangential to the boundary. Thus there are no ingoing constraint modes and 
constraint preservation is straightforward. 

In the strongly well-posed harmonic system described in Sec. [Xj the harmonic conditions C" became surrogates 
for the Hamiltonian and momentum constraints. Because C° satisfies a homogeneous wave equation, there are four 
ingoing constraint modes. These could be eliminated by dissipative boundary conditions with homogeneous data. In 
Sec. [Xj homogeneous Dirichlet conditions on C" were chosen. This allowed the constraints to be enforced in terms 
of first differential Sommerfeld conditions on the metric. Homogeneous Neumann or Sommerfeld conditions on the 
constraints would also ensure constraint preservation but at the expense of a more complicated coupling with the 
evolution system for the metric. 

The worldtube constraints which arise in the gravitational version of the null-timelike IBVP discussed in Sec. I VI CI 
present an entirely different aspect. In that problem, boundary data on a worldtube T and initial data on an outgoing 
null hypersurface Mo determine the exterior spacetime by integration along the outgoing nullgeodesics. The worldtube 
constraints impose conditions on the integration constants. The Bondi-Sachs formalism [3|, |59| introduces coordinates 
x" = {u,r,x ) based upon a family of outgoing null hypersurfaces Afu, where u labels the null hypersurfaces, x 
are angular labels for the null rays and r is a surface area coordinate. The evolution system is composed of radial 
propagation equations along the outgoing null rays consisting of the hypersurface equations G" = C^V yU — 0, which 
only contain derivatives tangent to A/'u, and the evolution equations Gab — hdAsg^^GcD — 0. 

The components of Einstein's equations independent of the hypersurface and evolution equations are worldtube 
constraints (called supplementary conditions by Bondi and Sachs), 

g^^'GAB = (11.1) 

G'a = (11.2) 

Gl = 0. (11.3) 

When the hypersurface and evolution equations are satisfied, the contracted Bianchi identity 

VyG""^ = (11.4) 

implies that these equations need only be satisfied on the worldtube T. The identity for fi — r reduces to q^^Gab = 
so that (jll.ip becomes trivially satisfied. (Here it is necessary that the worldtube have nonvanishing expansion so 
that the areal radius r is a non-singular coordinate.) The identity for v = A then reduces to the radial ODE 

drir^G'-A) = 0, (11.5) 

so that G\ vanishes if it vanishes on T- When G^ = 0, the identity ioi v = u then reduces to 

drir^Gl) = 0, (11.6) 

so that G^ also vanishes if it vanishes on T. 

Thus the worldtube constraints reduce to (|11.2p and ( 111.31 ). which are equivalent to the condition that the Einstein 
tensor satisfy 

eG^N, = 0, (11.7) 

where ^'' is any vector field tangent to the worldtube, whose normal is A^^. These are the boundary analogue of the 
momentum constraints for the Cauchy problem. In Stewart's treatment of the 3-1-1 IBVP (see Sec. lVII A|) . it was the 
Cauchy momentum constraints which were enforced on the boundary. In the characteristic IBVP, it is the worldtube 
constraints (|11.7p which must be enforced. They form three components of the boundary constraints (|7.ip proposed 
by Frittelli and Gomez for the 3-1-1 problem. (See Sec. I VII Cl ) 

This worldtube constraints (|11.7p can be interpreted as flux conservation laws for the ^-momentum contained in 
the worldtube fill. 



p^{u2) - Pd^i) = / dSA^'v^e - v^yf'-e'^)} (11. 

J Ui 



where 

P^= idS^M^C""^ (11.9) 
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and dSf^i, and dSi, are the 2-surface and S-volume elements on the worldtube. When £^^ is a Kilhng vector for the 
intrinsic 3-metric of the world-tube, this gives rise to a strict conservation law. For the limiting case at X"*", these 
flux conservation laws g overn the energy-momentum, angular momentum and supermomentum corresponding to the 
asymptotic symmetries [143|. For an asymptotic time translation, they give rise to the Bondi's famous result Q 
relating the mass loss to the square of the news function. 
In terms of the intrinsic metric of the worldtube 

H^, = g^, - N^^N,, (11.10) 

its intrinsic covariant derivative Vp, and its extrinsic curvature 

K^, = HPVpN,, (11.11) 

the worldtube constraints (|11.7p can be rewritten as 

Hl^G^pNP = V,{Kl^ - SiKP) = 0. (11.12) 

These are the analogue of (|2.17p for the Cauchy problem but they now form a symmetric hyperbolic system because 
of the timelike nature of the worldtube. In t erms of a dyad (|10.15p adapted to the foliation of the worldtube, this 
gives rise to the Worldtube Theorem [l44| : 

Given Hab, Q'^Q Kab and K , the worldtube constraints constitute a well-posed initial-value problem which determines 
the remaining components of the extrinsic curvature Kab- 

The theorem constrains the integration constants for the nullcone-worldtube IBVP. Similarly, they constrain the 
boundary data for a 3 -I- 1 IBVP subject to the Frittclli- Gomez conditions. Unfortunately, for neither of these IBVPs 
has it been possible to combine the boundary constraints with the evolution system in a manner consistent with a 
strongly well-posed IBVP. 

The enforcement of the boundary constraints is an indirect way to enforce the Hamiltonian and momentum con- 
straints constraints H = G'^'^n^n^ and P^ = h'^G'"^ni,, where in the 3 + 1 decomposition with respect to the Cauchy 
hypersurfaces 

The more direct approach commonly used to investigate constraint preservation in the 3 + 1 Cauchy problem is to cast 
the contracted Bianchi identity (|11.4p into a hyperbolic system. The results depend upon the particular formulation. 
As a first example, consider the ADM system (j2.15p in which only the 6 Einstein equations 

hP.h'^.R/ = (11.13) 

are evolved. Application of the contracted Bianchi identity gives rise to the symmetric hyperbolic constraint propa- 
gation system 

n'^d^P'-h'^djH = B^'^G^^n", (11.14) 

where the coefficients B'^ and -B^'*' arise from Christoffel symbols and do not enter the principal part. When applied 
to the IBVP, a complication arises from the component of the shift normal to the boundary, 

pN ^ ^^jy^ ^ _^ gjjjjj Q (11.15) 

in terms of the lapse a and the hyperbolic angle 9 p.4p governing the velocity of the boundary. Here f3^ < {f3^ > 0) 
for a boundary which is moving inward (outward) with respect to the Cauchy hypersurfaces. An analysis of (111.14^ 
shows that only one boundary condition is allowed provided /3^ < 0, i.e provided the boundary is moving inward. In 
that case, the theory of symmetric hyperbolic systems guarantees that all the constraints would be preserved if the 
single constraint 

H + P'N, = Gf^^n'^K" = (11.16) 

is satisfied at the boundary, where K'^ is the outgoing null vector to the foliation of the boundary. (Additional 
boundary conditions are necessary for constraint preservation if f3^ > 0.) By virtue of the evolution system (jll.131) . 
the constraint (|11.16l) is equivalent to 

Gni.Kf'K'' = 0. (11-17) 
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This is the Raychaudhuri equation (of. [l45J ) 

K^'d^,e + ^e^ + aa^O, (11.18) 

where 9 is the expansion and a is the shear of the outgoing null rays tangent to K^. Thus, for the ADM system, 
constraint preservation can be enforced by the Somnierfeld boundary condition (|11.18|) for 6. Unfortunately, although 
the constraint system has these attractive properties, the ADM evolution system is only weakly hyperbolic and 
consequently leads to unstable evolution. 

Next consider the BSSN evolution system, which enforces the 6 Einstein equations 



KKRp'' ~ -h;H = 0. (11.19) 



The contracted Bianchi identity now implies the constraint system 



rfd-^H-djP^ = B-fG^^n" (11.20) 

1 

— / 
3 



n^d^P' + -h^^djH = B'^'^G^^n'^. (11-21) 



This is not symmetric hyperbolic and would not lead to stable constraint preservation even for the Cauchy problem. 
This is remedied in the course of introducing auxiliary variables which reduce the BSSN system to first order form. 
Auxiliary constraints are mixed into the evolution system (J11.19P and they combine with the constraint system (|11.2ip 
to form a larger symmetric hyperbolic constraint system. There is a large freedom in the constraint-mixing parameters 
and gauge conditions. For a particular choice made by Nunez and Sarbach [75|, the linearization off Minkowski space 
yields a symmetric hyperbolic evolution system. The boundary conditions for this system are complicated by the 
normal component of the shift. As discussed in conjunction with p.Sp . the number of boundary conditions required 
by the advection equations introduced in the first order reduction depends upon whether (3^ is positive or negative. 
This forces use of a Dirichlet condition, e.g. /3^ = 0, rather than a Sommerfeld condition on the shift. Constraint 
preservation holds only in a certain parameter range, (6i < 1,62 < 1) for t he b oundary conditions given in equation 
(97) of [7a|. The particular choice 61 = 0, leads to the boundary condition [l46{ 

H - 3P'N, ^ Gi^^n^'in" - SN") = Z, (11.22) 

where Z represents contributions from the auxiliary constraints, or, by using the evolution system (I11.19p . 

G^,L^'L'' = Z, (11.23) 

where L'^ is the ingoing null vector to the boundary. It is a bizarre feature of the 3 + 1 problem that the constraint 
preserving boundary condition switches from the outgoing Raychaudhuri form (jll.lTp to the ingoing Raychaudhuri 
form (|11.23|) in going from the ADM to the BSSN system. The Raychaudhuri equation for the outgoing null direction 
cannot be imposed in the allowed range of (61, 62). 

The widely varying nature of constraint enforcement among different formulations does not provide any apparent 
insight. However, one problem common to many first order formulations arises from the advective derivative ri^S^, 
which determines whether the auxiliary variables are ingoing or outgoing at the boundary, depending on the sign of 
P^ . This problem could be avoided by instead using the derivative t^9^ determined by the evolution field, which 
can always be chosen tangential to the boundary. That suggests that the projection operator n^ associated with f, 
given in (14. 4p . might be useful in separating out the evolution system from the constraints, rather than the projection 
operator /i(^ used in (|11.13|) and (jll.l9|) . This gives rise to many ways to obtain a symmetric hyperbolic constraint 
system whose boundary treatment is independent of (3^ . As a simple example, in adapted coordinates the evolution 
system G*-' = A(5*^G**, with A > 0, leads via (|11.4p to the symmetrizable constraint system 

dtG*^ + djG*^ = lower order terms 
dtG*' + Aa,;G" = lower order terms. (11.24) 

Independently of /3^, this system requires only 1 boundary condition to preserve all the constraints. 
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XII. GEOMETRIC UNIQUENESS OF THE IBVP 

The solution of the Cauchy problem has the important property of geometric uniqueness, i.e. Cauchy data {hab, kab) 
on So determine a metric gab which is unique up to difFcomorphism. Under a diffeomorphism ip, the data {'4'*hab, '4'*kab) 
determines an equivalent metric. As well as being a pretty result, this has the practical application of allowing 
numerical simulations with the same initial data but carried out with different formulations and different gauge 
conditions to produce geometrically equivalent spacetimes. Friedrich [35] has emphasized that this property remains 
an unresolved issue for the IBVP. 

There are different ways in which this property might be formulated for the IBVP. The most demanding way 
would be to require that the data at a point of the boundary be locally determined by the boundary geometry in 
the neighborhood of that point. Such data might include the trace K of the extrinsic curvature of the boundary, 
which forms part of the data for the Friedrich-Nagy system. However, it is clear that at least two more pieces of 
data are necessary to prescribe the gravitational radiation degrees of freedom. In the Friedrich-Nagy system, these 
two pieces of data are supplied by the combination (|9.3p of the Weyl tensor components ^o and ^4. However, the 
associated outgoing and ingoing null vectors are not determined by the local geometry but depend upon the choice of 
timelike evolution field T"^ tangent to the boundary, according to (j9.4p . This could be avoided by requiring the se nu ll 
vectors to satisfy the local geometric condition that they be principle null directions of the Weyl tensor (cf. 145] ) ; 
but in a general spacetime this would lead to four choices which would then have to be incorporated somehow into 
the evolution system. An alternative, suggested in [35| , is to base the data on the eigenvector problem determined by 
the trace free part of the extrinsic curvature of the boundary, 

{Kab - \HabK)V'' = \HabV\ (12.1) 

As in the preceding presentation, Hab is the intrinsic metric of the boundary. For a spherical worldtube r — R m. 
Minkowski space, 

Kab - \HabK = -^{Hab + ifafb) (12.2) 

where Ta is a timelike eigenvector. This raises the possibility of whether this eigenvector problem can be used to pick 
out a locally preferred timelike direction Ta in the curved space case. Similar algebraic properties of the extrinsic 
curvature hold under roundness conditions which are typically satisfied by the artificial outer boundary of an isolated 
system. However, whether such an approach can be incorporated into the evolution system and whether the two 
radiation degrees of freedom can be encoded in the extrinsic curvature are not obvious. 

Neither of the two strongly well-posed formulations of the IBVP described in See's IIXI and 1x1 are based upon purely 
local geometric data. In both of them, a foliation of the boundary consistent with a choice of evolution field plays an 
essential nonlocal role. This suggests that a version of geometric uniqueness based upon purely local data might not 
be possible. The prescription of an evolution field i" as part of the boundary data provides the necessary structure 
to pose a version of geometric uniqueness JSJ, 13^. As explained in See's lIVl and FVl the fiow of the evolution field 
carries the initial edge Bq into a foliation Bt of the boundary; and it carries the initial Cauchy data into a stationary 
background metric ^^^ according to (J5.12I) . Thus the evolution field provides the two essential structures to geometrize 
the boundary data: the foliation Bt determines the outgoing null direction K"^ and the preferred background metric 
allows the Sommerfeld derivative to be expressed covariantly as K°'Va in terms of the background connection. Under 
a diffeomorphism, the evolution field transforms according to i° — >■ '0*^° with the consequence that g^^ -^ i'*§ab- 

The boundary data q° and a for the covariant version of the covariant Sommerfeld conditions p0.17p (|10.23p then 
have the geometric interpretation that 

q'' - K\^b - ^b)K'' (12.3) 

is the acceleration of the outgoing null vector K'^ relative to the background acceleration, and 

a = ig"Q''(V, - Va)Kb, (12.4) 

is the shear of K'^ relative to the background. The use of the shear in posing geometrical boundary conditions for 
the harmonic formulation was also suggested in (l02| . The rotation freedom in the dyad dependence (|12.4p can be 
removed by introducing the rank-2 shear tensor 

a-b ^ ^{Q-CQbd _ ^QabQcd^(^^^ _ vj^^ ^ ^ ^ QaQba"', (12.5) 
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with cr'^^Vbi = 0. 

By construction, all quantities involved in the boundary conditions map as tensor fields under a diffeoniorphism ip. 
As a result of the covariant form of the generalized harmonic equations ()10.10p and (jl0.11|) , the solution fab = gab — dab 
also maps as a tensor field. The metric gab satisfies the generalized harmonic condition (|10.12l) with respect to the 
background Qab and the mapped metric ii'*qn.b satisfies the generalized harmonic condition with respect to i^* gab 

This can be taken one step further [34| . As characterized in [4(| , the Cauchy data hab and kab can be interpreted 
as fields hab and kab on a disembodied 3-manifold Sq via its embedding Sq in the 4-dimensional spacetime manifold 
M. A similar approach applies to the boundary data. Let 

q" - qjvA^" + gf , (12.6) 

SO that q^ is tangent to the boundary T. Then the fields qN, Qj- and cr°^ are intrinsic to the boundary. Along with 
the Cauchy data and the hyperbolic angle Q at the edge Sq, they can be induced by the embedding of a disembodied 
version of data. This leads via the well-posedness of the IBVP with Sommerfeld data to a harmonic version of a 

Geometric Uniqueness Theorem: 

Consider the 3-manifolds T and So meeting in an edge Bq. On So prescribe the symmetric tensor fields hab and kab, 
subject to the Hamiltonian and momentum constraints and the condition that hab be a Riemannian metric . On Bq 
prescribe the scalar field 0. On T prescribe a smooth foliation Bt, parametrized by a scalar function t, the scalar 
field qN, the vector field qj- and the rank-2 tensor field ct"''. Then, after embedding 5o U T as the boundary SqUT 
of a 4-manifold A4 as depicted in Fig. HI this provides the Sommerfeld boundary data for a vacuum spacetime in a 
neighborhood of Bq which is unique up to diffeomorphism. 

Because the boundary data contain gauge information, this version of geometric uniqueness is weaker than for 
the Cauchy problem. It is an open question whether the boundary data can be prescribed purely in terms of local 
geometric objects j35| . Note that the data contains no information about the 3-metric of the boundary, not even that 
it is a timelike 3-manifold. The geometrical interpretation of the data involves the metric of the embedded spacetime, 
whose existence is in the content of the theorem. The diffeomorphism freedom lies in the freedom in the embedding 
and in the choice of evolution field i° , which determines the gauge and background geometry. See [37[ for a discussion 
of these issues in the context of linearized gravitational theory. 

The spatial locality of the solution can be extended, say to a boundary with spherical topology, by patching solutions 
together. However, the locality in time presents a more complicated problem regarding the maximal development of 
the solution. For instance, the time foliation might develop a gauge pathology which prematurely stops the evolution. 
That makes it unclear how the maximal development for the Cauchy problem, as constructed by Choquet-Bruhat 
[147| |. might be generalized to the IBVP. A restart of the evolution at an intermediate time in order to extend the 
solution would introduce a new gauge, new initial data, a new evolution field and thus a new background metric. A 
maximal development based upon the original background would have to be based upon a maximal choice of evolution 
field. ^_^ 

The geometric nature of the Sommerfeld conditions (|10.17p - (ll0.20p allows them to be formally applied to any 
metric formulation. In a 3 + 1 formulation, the metric has the decomposition 

5m"^ = -"m"'' + ^^N^ + Q^^, (12.7) 

where N^^ is the unit normal to the boundary which lies in the Cauchy hypersurfaces and, as before, Q^^ — QipQu) 
is the 2-metric intrinsic to its foliation Bt- The Sommerfeld boundary conditions (|10.17l) - (110. ISp and (|10.20l) supply 
boundary data for the N^N'^kf^^, Q^N'^k^^ and Q^Q^k^^ components of the extrinsic curvature of the Cauchy 
foliation. However, (110. 19p supplies the boundary data for the normal component of the shift, which for many 
3 + 1 formulations would require a Dirichlet condition that fixes its sign. The remaining Sommerfeld boundary 
conditions (|10.2ip - (|10.23l) . which enforce the harmonic constraints, would also require modification depending upon 
the particular 3 + 1 gauge conditions. See [T^ for a discussion relevant to the BSSN formulation. Numerical tests 
would be necessary to study whether application of (|10.17p - (|10.18p and (|10.20p would improve the performance over 
the present boundary treatment of 3 + 1 systems. 



XIII. NUMERICAL TESTS 



Post and Votta [25| have emphasized that "Verification and validation establish the credibility of code predictions. 
Therefore it's very important to have a written record of verification and validation results." The validation of a code 
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implies that its predictions are in accord with observed phenomena. For the present status of numerical relativity, 
in the absence of any empirical observations, the burden falls completely on verification. Post and Votta list five 
verification techniques: 

1. "Comparing code results with an exact answer". 

2. "Establishing that the convergence rate of the truncation error with changing grid spacing is consistent with 
expectations" . 

3. "Comparing calculated with expected results for a problem especially manufactured to test the code" . 

4. "Monitoring conserved quantities and parameters, preservation of symmetry properties and other easily pre- 
dictable outcomes" . 

5. "Benchmarking - that is, comparing results from those with existing codes that can calculate similar problems" . 

The importance of the first four techniques has now been recognized by most numerical relativity groups and their 
implementation in practice has improved the integrity of the field. Individual groups cannot easily carry out th e fifth 
technique independently. This was the motivation behind formation of the Apples with Apples ( AwA) Alliance [148| . 

The early attempts at developing numerical codes were primarily judged by their ability to simulate black holes, 
understandably because of the astrophysical importance of quantifying that system. When the difficulty with numer- 
ical stability became apparent, there was increased focus on a better mathematical and computational understanding 
of the analytic and numerical algorithms. Only a few groups had based their codes upon symmetric or strongly 
hyperbolic formulations of Einstein's equations and fewer had even begun to worry about how to apply boundary 
conditions. The cross fertilization between computational mathematics and numerical relativity was entering a pro- 
ductive stage. At the same time, standardized tests were developed by the AwA Alliance in order to isolate problems, 
calibrate accuracy and compare code results, http://www.ApplesWithApples.org, 

Such testbeds have been historically used in computational hydrodynamics. There are two fundamentally different 
types. One compares simulations of a physically important process, such as the binary black hole problem. The 
second type involve idealized situations which isolate problems, such as the "shock tube" test in computational fluid 
dynamics. This is the type of testbed considered by the AwA Alliance. 

The first tests were designed to study evolution algorithms in the absence of boundaries jl49l . Il50l | . Five tests were 
based upon a toroidal 3-manifold (equivalent to periodic boundary conditions): 

• The robust stability test evolves random initial data in the linearized regime. This is a pass/fail test designed 
as a screen to eliminate unstable codes. 

• The linearized wave test propagates a periodic plane wave either parallel or diagonal to an axis of the 3-torus. 
The test checks the accuracy in tracking both the amplitude and phase of the wave. 

• The gauge wave test is a pure gauge version of the linearized wave test, but with amplitude in the non-linear 
regime. 

• The shifted gauge wave test is based upon a gauge wave with non-vanishing shift. Both the gauge wave and 
shifted gauge wave tests are challenging because of exponentially growing modes in the analytic problem [3l[ . 

• The Gowdy wave test simulates an expanding or contracting toroidal spacetime, which contains a plane polarized 
gravitational wave in a genuinely curved, strong field context. 

The wave tests provide exact solutions which allow convergence measurements. Instabilities are monitored by the 
growth of the Hamiltonian constraint. Test results were carried out for code s ba sed upon numerous formulations: 
harmonic, Priedrich-Nagy, NOR, and several versions of BSSN and ADM. See |l5Cll | for the test results. 

Tests of the Cauchy evolution algorithm cull out algorithms whose boundary stability is doomed from the outset. 
A subsequent plan for boundary tests was formulated by opening up one axis of the 3-torus to form a manifold with 
boundary. This has the advantage that the boundaries are smooth 2-tori, thus avoiding the complication of sharp 
boundary points. This could later be extended to opening up all three axes to test performance with a cubic boundary. 
For the robust stability test the boundary data consist of random numbers. For the wave tests, the boundary data is 
supplied by the exact (or linearized) solution. See |l48J for the detailed specifications of the five AwA boundary tests. 

These boundary tests were first formulated and applied in the early development of boundary algorithms for 
harmonic codes. The robust stability test [7g was used to verify the stability of a code based upon a well posed 
IBVP for linearized harmonic gravity |77| . Subsequently, the gauge wave tests were carried out w ith a harmonic code 
whose underlying IBVP was well-posed for homogeneous Dirichlet and Neumann conditions |13l . Il03| . This revealed 
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problems in the very nonlinear regime, where the approximation of small boundary data was violated. The shifted 
gauge wave test po sed a n additional difficulty, beyond the unstable analytic modes that already challenged the Cauchy 
evolution test [3l|, Il5l| . The periodic time variation of the shift produced an effective oscillation of the boundaries 
which blue shifted and trapped the error resulting from the reflecting boundary conditions. This led to unstable 
behavior in the nonlinear regime. 

After formulation of the strongly well-posed harmonic IB VP descr ibed in Sec. |Xl the Sommerfeld boundary con- 
ditions were implemented and tested in a harmonic code |152h155| |. Numerical stability and accurate phase and 
amplitude tracking were confirmed by the robust stability and linearized wave tests. An important attribute of strong 
well-posedness is the estimate of boundary values provided by the energy conservation obeyed by the principal part of 
the equations. This boundary stability extends to the semi-discrete system obtained by replacing spatial derivatives 
by finite differences obeying summation by parts (the disc rete counterpart of integration by parts), so that energy 
conservation caries over to the semi-discrete problem [156j. Stability then ext ends to the fully discretized evolution 
algorithm obtained with an appropriate time integrator, such as Runge-Kutta [l57[. It was found that these discrete 
conservation laws were both effective and essential in controlling the exponential analytic modes latent in the gauge 
wave test. Although the analytic proof of well-posedness given in [18!| was based upon a scalar wave energy differing 
by a small boost from the standard energy, it is interesting that these successful code tests were based upon the 
standard energy. This confirms the robustness of the underlying approach. 

The oscillating boundaries in the shifted gauge wave test excite a different type of long wavelength instability which 
could not be suppressed by purely numerical techniques. Knowledge of the Sommerfeld boundary data allows the wave 
to enter and leave the boundaries, but the numerical error, although small and convergent, excites an exponential 
mode of the analytic problem. However, because this instability viola tes t he harmonic constraints it was possible to 
suppress it by a harmonic constraint adjustment of the form (j2.12p i3l|,|l52|. This example emphasizes the importance 
of understanding instabilities in the analytic problem in order to control them in a numeri cal s imulation. 

Other wave solutions which have been used for numerical tests are the Teukolsky waves llSSll . which are linearized 
spherical waves appropriate for testing a spherical boundary, and the nonlinear Brill waves [IGOJ ] , which are useful for 
testing wave propagation during collapse to a black hole. Teukolsky wave tests of the harmonic Sommerfeld conditions 
confirm that the constraint violation erro r due to homogeneous outer boundary data drops to numerical truncation 
error as the wave propagates off the grid |154| . Furthermore, when constraint damping is applied to the interior of 
the grid, the error drops to machine round-off. For Brill wave tests with the same code, homogeneous Sommerfeld 
conditions lead to considerable back reflection off the boundary, as expected from the discussion in Sec. IVIIII In 159] , 
Teukolsky waves were also used to te st an implementation of the improved higher order harmonic boundary conditions 
proposed by Buchman and Sarbach |lOd |. 

Rinne, Lindblom and Scheel [98| have developed a numerical test for comparing back reflection of waves from a 
spherical boundary. First, using perturbative techniques, they construct a reference solution for a linearized gravita- 
tional wave propagating on an exterior Schwarzschild background. Then the full numerical code is run with a finite 
spherical outer boundary. The error with respect to the reference sol ution is used to compare different choices of 
boundary conditions. Using the first order harmonic code described in [161], they used this test to compare several 
boundary conditions. The comparisons of homogeneous Sommerfeld boundary conditions were in accord with the 
theoretical expectations discussed in See's IVlIIl and 1x1 The higher order Sommerfeld conditions and the ^I'o freezing 
condition both produced less back refection than the first order conditions (|10.27l) - (|10.30l) . The test was also used 
to reveal the spurious effects arising from a sponge boundary condition and also from the combination of nume rical 



dissipation and spatial compactification used by Pretorius [28|, |29j . In an independent follow up of this test |102| , the 
advantages of higher order Sommerfeld conditions were confirmed. 

Although 3-1-1 codes have been predominant in binary black hole simulations, boundary tests have not been very 
extensive as compared to harmonic codes. This is especially pertinent for the BSSN system. Boundary tests based 
on a gravitational wave perturbation of the Schwarzschild exterior have been carried out for the KST system [96|. 
The results revealed instabilities although they showed how improvements to the boundary conditions could reduce 
constraint violation. The robust stability and Brill wave tests have been applied to boundary conditions for a 
symmetric hyperbolic version of the Einstein- Christoffel system [93] . Although the Cauchy problem for this system 
is well-posed, both tests revealed instabilities due to the boundary algorithm. Again the tests were useful guides for 
understanding the problems with the boundary treatment. 

C. Bona and C. Bona-Casas [ 162| have made the first application of the Gowdy wave boundary test. They show how 
it can be applied to a symmetric hyperbolic version of the first order Z4 formalism. Here, as in many other studies, 
the Cauchy problem is well-posed but the strong well-posedness of the IBVP depends on the boundary condition. 
Improper boundary conditions can lead to instability and/or constraint violation. They first demonstrate that the 
test is effective at investigating methods for preserving the energy constraint for the Z4 system in a strong field 
environment. In a subsequent work [l63| . they apply both the Gowdy wave and robust stability boundary tests in an 
expanded study of constraint violation in the Z4 framework. The robust stability case is applied both with opening 
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up one axis of the 3-torus and with a fully cubic boundary. This allows testing SBP algorithms at the corners and 
edges. The results show numerical stability of the proposed boundary algorithms in the linearized regime. The Gowdy 
wave test extends this study to constraint preservation in the nonlinear regime. The results provide further evidence 
of numerical stability and show that constraint violation can be kept at the level of discretization error. 

XIV. OPEN QUESTIONS 



The IBVP has analytical, computational, geometrical and physical aspects. The analytic goal is a strongly well- 
posed IBVP, which is the prime necessity for the computational goal of an accurate evolution algorithm. It is also the 
raison d'etre for the geometric goal of a gauge invariant formulation of the boundary data. The prime physical goal, 
at present, is the accurate simulation of binary black holes. The binary black hole problem has taken a course of its 
own, which has been remarkably successful in view of the gaps in our current understanding of the other aspects of 
the IBVP. 

Some important open questions which would help close those gaps are: 

• Question 1. Is there a strongly well-posed IBVP based upon a 3 + 1 formulation? 

Some insight into this question would be provided by the answer to 

• Question 2. Can the necessary boundary data be represented by gauge invariant, local geometric objects? 

In the Friedrich-Nagy treatment, there are three pieces of boundary data which are not pure gauge: the trace 
K of the extrinsic boundary curvature, which determines the location of the boundary, and the Weyl curvature 
components encoding the two radiation degrees of freedom. In the harmonic system, the Sommerfeld data which 
encode the radiation consist of the shear, or the Weyl curvature components for a higher order condition. This 
leads to 

• Question 3. In the harmonic formulation, can the trace K be used as gauge invariant boundary data? 
For the Cauchy problem, there exists a maximal development of the solution [l47|. 

• Question 4. What is the proper formulation of the maximal development of the IBVP? 

The quote of Turing at the beginning of this review expresses the relative degree of difficulty between the Cauchy 
problem and the IBVP. 
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